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I.  INTRODUCTION 


This  report  details  the  progress  made  by  ARAP  over  the  past  year 
towards  the  development  of  a  computer  model  for  determining  the 
detailed  low-level  atmospheric  distributions  of  velocity,  temperature, 
moisture,  refractive  index,  and  the  turbulent  variances  of  these 
quantities  for  marine  environments.  Since  low-level  clouds  or  fog  are 
a  frequent  occurrence  in  the  marine  boundary  layer,  the  prediction  of  the 
formation,  evolution,  and  dissipation  of  these  features  is  a  necessary 
integral  part  of  the  desired  model.  In  addition  to  appropriately 
modeling  the  turbulent  transport  of  momentum,  heat,  and  moisture,  this 
necessitates  the  incorporation  of  the  physics  of  thermal  radiation  and 
moisture  change  of  phase  into  the  boundary  layer  dynamics. 

In  the  next  section  we  give  a  very  brief  review  of  our  research 
program  in  marine  boundary  layer  modeling.  This  is  followed  by  a 
discussion  of  model  developments  and  model  calculations  performed  over 
the  period  of  this  contract.  Three  papers  prepared  for  publication 
during  this  period  are  incorporated  as  Appendices  A,  B  and  C. 


II.  BACKGROUND 


ARAP's  approach  to  the  problem  of  modeling  the  atmospheric 
boundary  layer  for  marine  environments  has  been  to  use  the  invariant 
second-order  closure  model  of  turbulence  developed  by  Dr.  Coleman  duP. 
Donaldson  and  his  associates  at  ARAP  over  the  past  few  years.  The 
fundamentals  of  this  approach  are  given  in  Reference  1 .  A  review  of 
the  status  of  this  model  as  applied  to  a  wide  variety  of  turbulent 
flow  problems  is  given  in  Reference  2.  Particular  applications  of  the 
model  as  applied  to  atmospheric  problems,  including  comparisons  with 
experimental  data,  are  documented  in  References  2-6. 

Reference  7  is  a  technical  report  detailing  the  model  development, 
sample  calculations,  and  verification  comparisons  made  under  our 
initial  contract.  It  describes  the  addition  of  humidity  and  the 
second-order  turbulence  correlations  involving  humidity  as  variables 
to  our  dry  atmospheric  boundary  layer  program  previously  developed  for 
the  Environmental  Protection  Agency.  Using  the  predicted  distributions 
of  temperature,  humidity,  and  pressure,  a  calculation  of  the  modified 
refractive  index,  M,  was  incorporated  in  the  program.  Local  minlmums 
In  the  M  distribution  with  respect  to  altitude  directly  indicate  the 
presence  of  a  radar  duct.  Since  we  are  predicting  the  second-order 
correlations  between  the  turbulent  fluctuations  in  temperature  and 
humidity  as  well  as  the  average  scale  of  the  turbulence,  we  have 
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available  the  information  to  also  compute  the  structure  of  the 
fluctuations  in  refractivity. 


Reference  7  includes  the  results  of  several  sample  calculations; 
e.g.:  (a)  a  sample  calculation  using  output  from  FNWC  (supplied  by 

J.  Kaitala)  as  upper  and  surface  boundary  conditions  for  our  boundary 
layer  program;  (b)  a  comparison  of  predicted  temperature  structure 
parameters  with  the  observations  of  Wyngaard,  et  al.  (Reference  8); 

(c)  a  calculation  with  boundary  conditions  roughly  corresponding  to 
the  conditions  observed  in  the  Atlantic  Tradewind  Experiment  by 
Augstein,  Schmidt,  and  Ostapoff  (Reference  9);  and  (d)  a  calculation 
simulating  shoreline  conditions  for  either  a  dry  land  breeze  over  the 
sea  or  a  moist  sea  breeze  over  dry  land. 

The  second  year  of  the  research  called  for  two  major  modifications 
to  the  model  described  in  Reference  7.  These  were  to  (a)  increase  the 
dimensions  of  the  program  to  a  two-dimensional,  unsteady  calculation 
to  permit  the  prediction  of  shoreline  conditions  developing  in  time, 
and  (b)  to  incorporate  the  radiative  flux  divergence  term  into  the 
one-dimensional  system  of  equations  in  a  coupled  manner.  These 
developments  are  detailed  in  Reference  10. 

The  two-dimensional,  unsteady  version  of  the  model  was  used  to 
calculate  the  typical  variation  in  the  coastal  planetary  boundary 
layer  (Reference  11).  The  resulting  diurnal  variation  in  the  sea 
breeze  induced  by  the  strong  stability  difference  in  the  boundary 
layer  response  over  the  land  and  that  over  the  water  produces  a  strong 
asymmetry  between  the  sea-breeze  and  the  land-breeze  circulation 
patterns.  In  previous  sea-breeze  models,  it  was  necessary  to  impose 
eddy  viscosities  which  were  a  strong  function  of  time  and  space  to 
gain  this  asymmetry.  In  the  present  model,  it  was  achieved  without 
the  need  to  introduce  any  new  empirical  information. 

The  incorporation  of  a  coupled  radiation  model  is  important  under 
stable  atmospheric  conditions  when  the  comparative  ratio  of  the 
divergence  of  radiation  heat  flux  to  that  of  the  turbulent  heat  flux 
may  reach  order  one.  The  primary  coupling  between  the  turbulent 
transport  and  radiation  comes  through  the  humidity  distribution.  The 
water  vapor  content  has  a  strong  influence  on  the  long-wave  radiative 
cooling,  while  the  liquid  water  content  strongly  controls  the  short¬ 
wave  radiative  heating.  Reference  10  describes  the  initial  radiative 
model  coupled  into  our  program  to  simulate  these  effects  in  the 
boundary  layer. 

The  third  year's  effort  was  divided  between  exercising  our  model 
for  verification  purposes,  and  extending  its  capability.  The  sample 
calculations  are  detailed  In  References  12  and  13.  These  sample 
calculations  demonstrate  the  strong  dynamic  interaction  between 
turbulent  transport  and  thermal  radiation.  In  general,  the  agreement 
between  model  predictions  and  field  observations  is  encouraging. 
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A  number  of  extensions  to  the  model's  capability  were  also  made. 
These  include:  (1)  the  removal  of  the  hydrostatic  approximation; 

(2)  the  incorporation  of  condensation  and  radiation  into  the 
two-dimensional,  unsteady  model;  and  (3)  the  ability  to  determine  the 
influence  of  particle  size  on  the  turbulent  transport  of  aerosols. 

We  also  investigated,  analytically,  some  of  the  relationships  between 
parameters  governing  distinct  fronts  in  the  planetary  boundary  layer, 
as  an  aid  toward  the  detailed  computation  of  the  structure  of  such 
fronts . 

Probably  the  most  important  result  detailed  in  the  most  recent 
contract  report  (Reference  14)  was  the  calculation  of  the  longitudinal 
roll-vortex  structure  which  frequently  characterizes  the  large 
turbulent  eddies  in  the  marine,  atmospheric  boundary  layer.  Results 
of  this  calculation  are  also  detailed  in  References  15  and  16.  From 
these  results  it  appears  that  in  the  middle  region  of  the  mixed  layer, 
more  than  half  of  the  vertical  momentum,  heat,  and  humidity  flux  is 
carried  by  the  large  roll  structure  rather  than  the  small-scale 
turbulence. 

Limited  comparisons  with  AMTEX  and  Great  Lakes  data  are  included 
in  Reference  14  along  with  a  detailed  discussion  of  the  field 
observation  requirements  for  definitive  model  verification  and  a 
discussion  of  several  improvements  in  the  numerical  techniques  used. 

Appendices  A,  B  and  C  are  three  papers  written  during  the 
present  contract  period.  Appendix  A  is  the  result  of  the  major  part 
of  our  effort  devoted  to  the  incorporation  of  an  anisotropic  length 
scale  into  our  boundary  layer  model.  It  is  being  submitted  to  the 
Journal  of  Fluid  Mechanics.  Appendix  B  was  written  for  a  workshop  on 
water  vapor  in  the  atmosphere  held  at  Vail,  Colorado,  in  September 
1979.  It  reviews  the  problem  of  modeling  the  transport  of  water 
vapor  through  the  atmospheric  boundary  layer.  It  will  be  published 
by  Academic  Press  later  this  year.  Appendix  C  is  being  published  as 
a  paper  in  a  volume  on  Turbulent  Shear  Flows  resulting  from  the  2nd 
International  Symposium  on  that  subject  held  in  London,  July  1979. 

This  is  an  expanded  version  of  the  paper  which  appeared  in  the 
Symposium  proceedings  and  which  was  included  as  an  appendix  to 
Reference  14. 


III.  MODEL  DEVELOPMENTS 


A.  Incorporation  of  Two  Scales 

A  major  part  of  our  effort  during  the  current  contract  has  been 
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devoted  to  the  incorporation  of  an  anisotropic  length  scale  into  our 
boundary  layer  model.  Physically  the  modification  amounts  to 
recognizing  that  when  there  is  a  significant  difference  between  the 
horizontal  scale  of  a  turbulent  fluctuation  and  its  vertical  scale 
two  modifications  in  the  model  are  indicated.  First,  the  pressure 
scrambling  term  which  we  have  previously  modeled  as  a  simple  tendency- 
towards-isotropy  term  must  be  broken  into  two  terms.  A  term  which 
tends  to  transfer  energy  from  the  component  with  the  shorter  length 
scale  to  that  with  the  longer  length  scale  must  be  added  to  the  simple 
tendency-towards-isotropy  term.  Second,  the  dissipation  length  scale 
needs  to  be  a  weighted  average  between  the  two  scales.  The  weighting 
we  have  chosen  is  the  inverse  of  the  cube  of  the  respective  velocity 
variances  since  the  dissipation  has  dimensions  of  a  velocity  cubed 
over  the  length  scale.  As  detailed  in  Appendix  A,  the  modification 
in  the  model  makes  a  major  improvement  in  predicting  the  horizontal 
velocity  variance  near  the  surface  under  unstable  conditions.  Our 
standard  model  leaves  aH/u*  as  a  function  of  Z/L  in  the  surface 
layer  with  no  direct  dependence  on  Z-j/L  .  Observations,  and  the 
modified  model,  show  that  just  the  opposite  occurs. 


B.  Comments  on  an  Integral  Model  Based  on  the  Second-Order  Closure 
Equations  ~  ~  ~  ~  ~ 

Considerable  time  was  spent  during  this  contract  period 
investigating  the  development  of  a  numerical  model  of  the  planetary 
boundary  layer  based  on  utilizing  integrals  of  the  second-order 
equations.  The  idea  is  to  make  maximum  use  of  a  few  vertical  grid 
points.  To  accomplish  this  it  is  necessary  to  position  the  grid 
where  the  variables  show  the  maximum  variation  and/or  make  use  of 
layer-averaged  equations  which  are  less  sensitive  to  local  gradients. 
We  have  attempted  two  different  approaches  to  this  problem.  One  is 
to  introduce  parameterized  vertical  distributions  for  the  mean 
variables  and  turbulent  correlations  of  interest  and  formally 
integrate  the  governing  equations,  including  the  Reynolds  stress 
equations  as  well  as  conservation  of  mean  momentum,  energy  and 
species.  Preliminary  results  of  this  approach  were  given  in  our  last 
status  report  (Reference  14).  It  shows  promise,  but  has  the 
inherent  problem  that  results  are  sensitive  to  the  parameterization 
chosen  and  when  appropriately  general  vertical  distributions  are 
chosen  the  resulting  equations  become  algebraically  cumbersome.  In 
essence  this  approach  reduces  the  work  of  the  computer  but  at  the 
expense  of  the  work  required  for  the  programmer.  In  the  past  year, 
we  have  attempted  to  modify  the  approach  such  that  it  shifts  more 
of  the  work  load  back  to  the  machine.  The  governing  equations 
were  transformed  to  a  variable  grid  coordinate  system  so  that 
the  grid  points  utilized  would  remain  at  fixed  fractions  of  the 
boundary  layer  as  the  boundary  layer  thickness  evolves  in  time  and 
space.  Thus,  even  under  meteorological  conditions  that  lead  to 
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vastly  different  boundary  layer  heights  all  of  the  grid  points  will 
remain  useful,  helping  to  define  the  boundary-layer  flow.  Integral 
conditions  are  used  to  determine  the  evolution  of  the  boundary  layer 
thickness . 

Numerical  tests  have  shown  that  five  grid  points  are  sufficient 
to  quite  reasonably  represent  a  neutral  planetary  boundary  layer  when 
the  correct  boundary-layer  thickness  is  prescribed  and  the  bottom 
boundary  condition  is  applied  at  z  =  0.15  where  the  boundary 
conditions  impose  a  matching  to  the  logarithmic  surface  layer.  Based 
on  this,  a  technique  has  been  formulated  where  a  surface  layer  is 
used  at  the  bottom  of  the  boundary  layer  and  an  inversion  layer,  if 
necessary,  added  at  the  top  so  that  the  finite  difference  equations 
are  only  required  to  resolve  the  central  bulk  of  the  boundary  layer. 
The  grid  points  in  this  finite  difference  region  are  to  be  located 
at  fixed  fractions  of  5  .  The  thickness  5  is  to  be  determined  by 
integral  constraints  taken  across  the  complete  boundary  layer. 

To  this  point  we  have  only  programmed  the  integral  constraints 
in  the  1-D,  neutral  case.  We  want  to  develop  a  completely  stable 
algorithm  for  this  limiting  case  before  attacking  the  general 
problem.  We  have  attempted  several  different  formulations  for  this. 
Currently,  the  most  promising  is  based  on  satisfying  the  integrated 
momentum  and  total  kinetic  energy  equations.  In  the  limiting  case 
these  integral  equations  may  be  written  as: 

^  =  -  u?  +  f V»  m 
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The  integrals  are  composed  of  a  contribution  from  the  surface  layer 
where  U,  V  and  q  are  given  by  the  asymptotic  form  derived  in 
Reference  17, 

U  =  ^  £n  z/z0 
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and  a  contribution  from  the  remainder  of  the  boundary  layer  where 
the  variables  are  defined  by  the  finite  difference  equations.  With 
the  evolution  of  the  integrals  determined  by  equations  (1)-  (3)  the 
boundary  layer  thickness  may  be  determined  from  the  equation 
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5  =  H (U  q  +  V^)/EA 


when  the  shape  parameter  H  is  known.  We  have  not  yet  found  a 
satisfactory  stable  algorithm  to  determine  H  while  letting  it  vary 
in  response  to  changing  conditions. 

Further  work  is  required  to  develop  a  stable  model  which  is 
generalizable  to  the  more  interesting  cases  with  temperature 
gradients  and  humidity,  but  we  remain  optimistic  that  the  combination 
of  integral  constraints  and  a  variable  grid  tied  to  the  boundary-layer 
thickness  will  yield  an  efficient  numerical  model. 


IV.  COMPARISONS  BETWEEN  THE  MODEL  PREDICTIONS 
OBTAINED  FROM  DR.  BURK'S  CURRENT  MODEL  AT 
NEPRF  AND  ARAP'S  ONE-DIMENSIONAL  MODEL 


Dr.  S.  Burk  has  a  one-dimensional,  atmospheric  boundary  layer 
model  running  at  NEPRF  which  uses  a  second-order  closure  turbulence 
model  similar  to  ARAP's  model,  but  differing  in  some  details.  An 
attempt  was  made  to  make  detailed  comparisons  between  the  predictions 
obtained  from  these  models.  Ideally  these  comparative  calculations 
should  have  been  carried  out  for  conditions  where  reliable  data  on  the 
"correct"  results  were  available.  This  was  done  only  to  a  quite 
limited  extent. 

The  first  comparison  was  for  the  case  of  free  convection  where 
comparison  with  laboratory  data  is  available.  Figures  (1)  and  (2) 
show  the  results  for  the  vertical  and  horizontal  velocity  variance. 

For  our  ARAP  model  both  the  2-scale  and  the  single-scale  results  are 
shown.  The  NEPRF  model  should  more  nearly  agree  with  the  single¬ 
scale  result.  The  nearly  20%  departure  of  NEPRF' s  result  for  the  maxi¬ 
mum  value  of  aw  from  that  predicted  by  ARAP's  model,  which  agrees  with 
the  data,  is  somewhat  surprising.  This  is  particularly  true  since 
Mellor  and  Yamada  (Reference  18)  have  previously  reported  results  from 
their  model,  upon  which  Burk's  model  is  based,  to  be  in  agreement  with 
the  data.  The  horizontal  variance  of  both  the  ARAP  single-scale  model 
and  the  NEPRF  model  suffer  the  same  deficiency  at  the  surface. 

Figure  (3)  compares  the  temperature  variance  for  the  different 
models  with  the  laboratory  measurements.  All  of  the  models  agree 
reasonably  with  the  data  in  the  bottom  two-thirds  of  the  layer.  In 
the  neighborhood  of  the  inversion  the  ARAP  two-scale  results  came 
closest  to  the  data  but  even  it  is  somewhat  lower  than  the  observations. 
In  this  region  a  large  part  of  the  temperature  fluctuation  is 
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2- scale  ARAP 


Figure  1  -  Vertical  velocity  variance  as  a  function  of  normalized 
height  in  a  free  convection  layer  as  predicted  by  three  models. 
(See  Appendix  A  for  a  more  complete  description  of  the  difference 
between  the  two  ARAP  models.) 
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Figure  2  -  Horizontal  velocity  variance  as  a  function  of  normalized 
height  in  a  free  convection  layer  as  predicted  by  three  models. 

(See  Appendix  A  for  a  more  complete  description  of  the  difference 
between  the  two  ARAP  models.) 


Figure  3  -  Temperature  variance  as  a  function  of  normalized  height  in 
a  free  convection  layer  as  predicted  by  three  models.  (See  Appendix  A 
for  a  more  complete  description  of  the  difference  between  the  two 
ARAP  models. ) 


associated  with  internal  gravity  waves  which  cannot  be  adequately 
represented  by  the  1-D,  turbulent  model  . 

The  other  cases  chosen  for  comparison  were  attempts  to  simulate 
observations  for  two  different  days  for  which  data  were  obtained 
during  the  multi-platform  Cooperative  Experiment  in  West  Coast 
Oceanography  and  Meteorology  in  1976  (CEWCOM-76)  performed  off  the 
coast  of  southern  California  from  26  September  to  12  October 
(Reference  19).  Adequate  data  was  not  available  to  make  a  rigorous 
comparison  between  model  predictions  and  observations  as  was  done 
in  the  free  convection  case.  Here  we  will  show  comparisons  between 
the  two  model  results  (the  ARAP  model  was  run  only  in  its  single¬ 
scale  version  in  this  case)  and  give  some  indication  of  how  sensitive 
the  results  were  to  boundary  conditions  which  were  not  completely 
specified  by  the  data  available  to  us. 

Figures  (4)  -  (15)  compare  the  results  of  the  two  model 
predictions  for  the  3-4  October  case.  Both  models  are  initialized 
with  the  profiles  available  from  the  1850  IT  sounding  [Figure  (16)]. 
Subsequently  the  sea  surface  time  variation  was  input  as  a  boundary 
condition  to  the  model  from  the  reported  data.  The  geostrophic  wind 
was  held  fixed  at  5  m/sec,  the  upper  level  lapse  rate  was  set  at 
2.5°K/Km  and  a  subsidence  velocity  was  set  to  vary  linearly  with 
height  (proportionality  factor  =  -  2  x  10'5  sec-1 ) .  Latitude  was  set 
at  32.75°  and  the  solar  angle  was  set  to  correspond  to  1850  LT  with 
a  declination  of  0°.  The  vertical  distribution  of  mean  variables  and 
primary  turbulent  variances  and  covariances  are  shown  for  three  times 
(0350,  1545  and  1850)  on  October  4th. 

The  two  models  give  generally  similar  results,  but  there  are 
noteable  quantitative  differences.  ARAP's  model  results  show  a 
somewhat  more  rapid  drop  in  mixed  layer  thickness  followed  by  little 
change  between  the  last  two  time  breaks.  The  NEPRF  results  show  a 
thicker  mixed  layer  at  the  first  time  break  and  smaller  thicknesses 
at  the  last  two  time  breaks.  The  height  of  the  inversion  as  predicted 
by  the  ARAP  model  is  a  somewhat  better  representation  of  the  reported 
inversion  height  variation,  but  the  mixed  layer  temperature  variation 
predicted  by  the  NEPRF  model  is  closer  to  the  observation.  More 
fundamentally,  both  models  predict  clouds  to  occur  as  given  by  the 
liquid  water  distribution  in  Figure  (15),  while  4  October  was  reported 
as  a  clear  day.  The  ARAP  model  predicts  denser  clouds.  Cloud 
evolution  is  somewhat  unstable  in  that  once  clouds  appear  they  lead 
to  increased  cooling  in  the  mixed  layer  through  cloud  top  radiation. 

This  increased  cooling  leads  to  increased  cloudiness. 

The  radiational  properties  of  the  upper  atmosphere  were  found  to 
be  quite  important  in  determining  the  results  of  this  boundary  layer 
run.  Figures  (17)  -  (22)  show  the  results  obtained  at  1545  when  the 
previous  run  was  repeated  with  only  one  change:  that  of  the  water 
vapor  content  in  the  total  atmosphere  above  the  top  of  the  computational 
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Figure  4  -  Comparison  between  model  predictions  for  mean  virtual 
potential  temperatures  at  three  different  times  on  3  -  4  October 
( -  ARAP  model,  .  NEPRF  model). 
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Figure  5  -  Comparison  between  model  predictions  for  mean  wind  in  the 
geostrophic  direction  at  three  different  times  on  3  -  4  October 
( -  ARAP  model,  .  NEPRF  model). 
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Figure  8  -  Comparison  between  model  predictions  for  total  velocity 

variance  at  three  different  times  on  3  -  4  October  ( -  A RAP  model, 

.  UEPRF  model). 
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Figure  9  -  Comparison  between  model  predictions  for  covariance  of  u 
and  w  velocities  at  three  different  times  on  3  -  4  October 
( -  ARAP  model, . NEPRF  model). 
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Figure  10  -  Comparison  between  model  predictions  for  covariance  of  v 
and  w  velocities  at  three  different  times  on  3  -  4  October 
( -  ARAP  model,  - NEPRF  model). 
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Figure  11  -  Comparison  between  model  predictions  for  covariance  of 
and  the  virtual  potential  temperature  at  three  different  times  on 
3-4  October  ( -  ARAP  model, . NEPRF  model). 


Figure  12  -  Comparison  between  model  predictions  for  covariance  of  w 
and  the  total  water  mixing  ratio  at  three  different  times  on 
3-4  October  ( -  ARAP  model, - '1EPRF  model). 
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Figure  13  -  Comparison  between  model  predictions  for  vertical  velocity 

variance  at  three  different  times  on  3  -  4  October  ( -  ARAP  model, 

.  NEPRF  model). 
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Figure  14  -  Comparison  between  model  predictions  for  variance  of  total 
water  mixing  ratio  at  three  different  times  on  3  -  4  October 
( -  ARAP  model,  .  NEPRF  model). 
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Figure  15  -  Comparison  between  model  predictions  for  mean  liquid 

water  content  at  three  different  times  on  3  -  4  October  ( -  ARAP 

model,  .  NEPRF  model). 
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Figure  16  -  Potential  temperature,  virtual  potential  and  water  vapor  mixing  ratio  used  to 
initialize  the  run  for  the  predictions  given  in  Figures  4-15. 


Figure  17  -  Comparison  between  model  predictions  for  mean  virtual 
potential  temperature  at  one  time  for  the  same  conditions  as  in 
Figure  4,  except  for  an  increase  in  upper  atmospheric  water  vapor. 
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Figure  18  -  Comparison  between  model  predictions  for  mean  wind  in  the 
geostrophic  direction  at  one  time  for  the  same  conditions  as  in 
Figure  5,  except  for  an  increase  in  upper  atmospheric  water  vapor. 
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Figure  19  -  Comparison  between  model  predictions  for  mean  wind 
perpendicular  to  the  geostrophic  direction  at  one  time  for  the  same 
conditions  as  in  Figure  6,  except  for  an  increase  in  upper  atmospheric 
water  vapor. 


Figure  20  -  Comparison  between  model  predictions  for  mean  total  water 
mixing  ratio  at  one  time  for  the  same  conditions  as  in  Figure  7, 
except  for  an  increase  in  upper  atmospheric  water  vapor. 


Figure  21  -  Comparison  between  model  predictions  for  total  velocity 
variance  at  one  time  for  the  same  conditions  as  in  Figure  8,  except 
for  an  increase  in  upper  atmospheric  water  vapor. 


Figure  22  -  Comparison  between  model  predictions  for  variance  of 
total  water  mixing  ratio  at  one  time  for  the  same  conditions  as  in 
Figure  14,  except  for  an  Increase  in  upper  atmospheric  water  vapor. 
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domain.  A  change  in  the  parameter  HU  (defined  as  the  average  mass 
density  of  H20  vapor  in  the  atmosphere  above  the  computational 
domain  divided  by  the  surface  density  of  H20)  from  10"3  to  10_1 
essentially  wipes  out  all  cloudiness.  It  also  reduces  the  inversion 
height  by  a  factor  of  2  to  3.  We  could  have  used  this  parameter  to 
try  to  make  the  model  results  match  the  limited  available  data  but 
this  was  not  the  intent.  We  are  interested  in  discerning  differences 
between  the  two  models. 

Figures  (23)  -  (34)  compare  the  results  for  9  October  with  the 
0350  PDT  sounding  used  to  initialize  the  run.  In  this  case  the 
subsidence  was  set  equal  to  zero  and  the  geostrophic  wind  was  held 
constant  at  8  m/sec  while  the  sea  surface  temperature  dropped  1.3°C 
over  a  14  hour  period.  The  upper  level  humidity  parameter  HU  was 
held  at  10“3.  The  two  model  predictions  agree  somewhat  better  in 
this  case  than  in  the  last.  Again  the  ARAP  model  leads  to  somewhat 
denser  clouds,  and  the  accompanying  effects  of  increased  cooling, 
increased  instability,  and  increased  turbulence.  At  1610  the 
temperature  variation  shown  for  the  ARAP  model  appears  to  be 
somewhat  closer  to  the  reported  data  with  a  mixed  layer  vertical 
potential  temperature  of  approximately  289  and  an  inversion  height  of 
approximately  280  m. 

In  summary  the  NEPRF  model  and  the  ARAP  1-D,  single-scale 
model  give  generally  similar  predictions.  However,  quite  significant 
quantitative  differences  do  occur,  particularly  in  cases  involving 
clouds.  In  all  the  cases  tested,  ARAP's  model  lead  to  higher  values 
of  liquid  water  content.  The  data  compared  with  herein  is  not 
adequate  to  clearly  distinguish  which  is  more  accurate. 


V.  RECOMMENDATIONS 


The  current  objectives  of  our  research  are  threefold.  They  are: 
(1)  to  improve  the  basic  physics  of  the  model;  (2)  to  exercise  the 
existing  model  to  exemplify  phenomena  which  may  be  predicted  to  occur 
and  to  test  model  validity;  and  (3)  to  simplify  the  numerics  to  make 
the  model  more  convenient  to  use.  Five  specific  tasks  have  been 
chosen  for  the  next  contract  period  to  provide  a  balanced  program 
toward  meeting  these  objectives. 

The  first  task  is  to  improve  the  cloud  physics  of  the  model.  The 
current  model  is  based  on  quite  simple  physics.  Thermodynamic 
equilibrium  is  assumed  to  exist  between  liquid  and  gas  phase  water  at 
all  times.  The  liquid  which  exists  is  assumed  to  be  in  the  form  of 
small  droplets  of  specified  size.  The  droplet  size  is  constant  in 
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Figure  23  -  Comparison  between  model  predictions  for  mean  virtual 
potential  temperatures  at  three  different  times  on  9  October 
( - -  ARAP  model, - NEPRF  model). 


Figure  24  -  Comparison  between  model  predictions  for  mean  wind  in 
the  geostrophic  direction  at  three  different  times  on  9  October 
r -  ARAP  model,  .  NEPRF  model). 
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Figure  25  -  Comparison  between  model  predictions  for  mean  wind 
perpendicular  to  the  geostrophic  direction  at  three  different  times 
on  9  October  { -  ARAP  model, - .NEPRF  model). 
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Figure  26  -  Comparison  between  model  predictions  for  mean  total 

water  mixing  ratic  at  three  different  times  on  9  October  ( -  ARAP 

model, - :1EPRF  model). 
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Figure  27  -  Comparison  between  model  predictions  for  total  velocity 

variance  at  three  different  times  on  9  October  ( -  ARAP  model, 

.  NEPRF  model). 
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Figure  28  -  Comparison  between  model  predictions  for  covariance  of  u 

and  w  velocities  at  three  different  times  on  9  October  ( -  ARAP 

model,  -  NEPRF  model). 
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Figure  31  -  Comparison  between  model  predictions  for  covariance  of  w 
and  the  total  water  mixing  ratio  at  three  different  times  on  9  October 
( -  ARAP  model, . MEPRF  model). 


Figure  32  -  Comparison  between  model  predictions  for  vertical  velocity 

variance  at  three  different  times  on  9  October  ( -  ARAP  model, 

-  MEPRF  model). 
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Figure  33  -  Comparison  between  model  predictions  for  variance  of 
total  water  mixing  ratio  at  three  different  times  on  9  October 
( -  ARAP  model, . NEPRF  model). 
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Figure  34  -  Comparison  between  model  predictions  for  mean  liquid 

water  content  at  three  different  times  on  9  October  ( -  ARAP  model, 

.  NEPRF  model). 
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space  and  time.  Only  the  number  density  varies  as  the  liquid  water 
content  varies.  In  actuality,  we  would  expect  the  droplet  size 
distribution  to  be  controlled  by  a  complex  interaction  between  the 
turbulent  fluctuations  in  relative  humidity  and  the  ambient  concentra¬ 
tion  of  condensation  nuclei.  Since  the  droplet  size  distribution  is 
such  an  important  variable  in  determining  visibility  within  a  fog  of 
given  liquid  water  content,  we  would  like  to  incorporate,  at  least, 
some  droplet  growth  dynamics  within  our  model.  We  will  attempt  to 
make  use  of  the  analytical  and  experimental  studies  performed  by 
other  NASC  contractors  on  fog  droplet  dynamics.  The  important 
interaction  between  thermal  radiation  and  droplet  size  will  also  be 
included. 

The  second  task  is  to  continue  development  of  the  integral  model. 
The  biggest  disadvantage  of  our  ARAP  model  is  its  numerical  complexity. 
It  requires  both  a  considerable  amount  of  computer  time  and  a 
considerable  degree  of  indoctrination  for  a  new  individual  to  make 
appropriate  choices  of  all  of  the  input  variables.  Our  current 
investigation  of  the  construction  of  a  simpler  numerical  model, 
based  on  utilizing  integrals  of  the  second-order  closure  equations, 
has  convinced  us  that  this  method  has  the  potential  to  make  our  model 
significantly  more  convenient  to  use  with  little  loss  in  accuracy. 

From  a  practical  point  of  view,  it  appears  that  the  accuracy  of  the 
simpler  model  will  generally  be  controlled  more  by  uncertainties  in 
the  meteorological  inputs  than  in  the  differences  in  model  errors 
between  the  simpler  model  and  the  full  model. 

The  third  task  is  to  make  comparisons  between  the  results  of 
other  mixed  layer  models,  such  as  Lilly's  (Reference  20),  and  the 
results  from  both  our  full  second-order  closure  model  and  our 
integral  model  to  be  developed  under  the  second  task.  Several 
parameterized  models  of  the  unstable  marine  boundary  layer  have 
appeared  in  the  literature.  Vie  propose  to  examine  the  differences 
predicted  by  the  most  promising  of  these  and  the  predictions  of  our 
full  second-order  closure  model.  This  will  permit  a  checking  of 
such  things  as  the  entrainment  parameterization  used  in  the  models. 

The  fourth  task  is  to  perform  specific  calculations  with  the 
existing  program  to  examine  the  different  characteristics  between 
warm-water  fog  and  cold-water  fog.  We  have  used  our  model  to  perform 
example  calculations  involving  fog  (Reference  12),  but  have  not  made 
a  systematic  investigation  of  the  differences  between  warm-water  and 
cold-water  fogs.  Observations  seem  to  suggest  that  warm-water  fogs 
are  found  more  often  than  cold-water  fogs.  Is  this  an  anomoly  of  the 
observations  or  is  there  a  fundamental  bias  in  the  physics  which  makes 
the  warm-water  fog  more  persistent?  By  using  our  model  to  examine  the 
interaction  between  radiative  and  turbulent  transport  within  a  fog 
bank,  we  should  be  able  to  answer  this  question. 

It  appears  appropriate  at  this  time,  after  more  than  five  years 
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of  development,  to  write  a  critical  review  of  our  ARAP  marine  boundary 
layer  model.  The  plan  under  Task  5  would  be  to  review  the  model 
formulation  and  verification  tests,  and  detail  the  current  model's 
capabilities  and  weaknesses.  This  task  would  not  lead  directly 
towards  the  accomplishment  of  any  of  the  three  previously  stated 
objectives  but  it  should  enhance  the  transferability  of  our  researcn 
results  to  others. 
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APPENDIX  A 

Incorporation  of  an  Anisotropic  Scale 
into  Second-Order  Closure  Modeling  of  the 
Reynolds  Stress  Equation 


by 

W.  S.  Lewellen  and  G.  Sandri 


(To  be  submitted  for  publication  to  the  Journal  of 
Fluid  Mechanics.) 


1 .  INTRODUCTION 


Modeling  turbulent  flow  is  a  formidable  challenge. 

Although  the  Navier-Stokes  equations  provide  a  true  representa¬ 
tion  of  an  incompressible  viscous  fluid,  the  wide  range  of 
scales  which  occur  at  high  Reynolds  numbers  serves  to  preclude 
using  these  equations  to  represent  the  flow  in  direct  calcula¬ 
tions  except  in  rare  instances.  A  variety  of  approaches  have 
been  offered  by  researchers  to  circumvent  this  problem  with 
various  degrees  of  success.  An  approach  which  has  yielded  at 
least  some  degree  of  success  has  been  to  extend  the  number  of 
variables  which  describe  the  fluid  to  include  the  second-order 
turbulent  correlations  of  the  primary  variables.  This  extended 
set  of  correlations  is  governed  by  an  extended  set  of  conserva¬ 
tion  laws,  the  Reynolds  stress  equations,  obtained  by  taking 
exact  second-order  moments  of  the  instantaneous  Navier-Stokes 
equations.  Second-order  closure  is  obtained  by  modeling  the 
higher-order  correlations  which  appear  in  the  Reynolds  stress 
equations  in  terms  of  the  describing  field  of  first-  and  second- 
order  correlations. 

We  will  not  attempt  to  give  a  review  of  all  of  the  models 
offered  in  the  literature  to  close  the  Reynolds  stress  equations. 
Past  reviews  have  been  given  by  Reynolds  (1976)  ,  Launder  and 
Spalding  (1972),  Lewellen  (1977),  Mellor  and  Herring  (1973). 
Although  no  complete,  up-to-date  review  of  this  research  area 
is  available,  we  believe  the  model  offered  herein  is  the  first 
attempt  to  incorporate  some  influence  of  an  anisotropic  scale 
into  a  numerical  model.  The  investigation  by  Hanjalic,  Launder 
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and  Schiestel  (1979)  of  multiple-t ime-scale  concepts  is 
probably  the  most  closely  related  work  appearing  in  the  liter¬ 
ature  . 

There  is  a  vast  amount  of  information  contained  in  the 
details  of  a  typical  turbulent  motion  because  of  the  wide  range 
of  scales  exhibited.  The  process  of  ensemble  averaging  neces¬ 
sarily  smooths  over  a  lot  of  this  information.  But  it  is  a 
prediction  of  the  means,  variances  and  covariances  of  the 
primary  variables  which  is  generally  most  desired  in  engineering 
models.  The  desire  is  to  provide  a  model  which  provides  an 
invariant  relationship  between  these  variables  for  turbulent 
flows  with  different  boundary  conditions.  It  appears  unlikely 
that  relationships  between  the  first-  and  second-order  moments 
at  a  single  point  can  be  uniquely  defined  without  resorting  to 
information  about  two-point  averages,  since  turbulence  is  a 
property  of  the  macroscopic  flow  field  rather  than  a  local 
property  of  the  fluid.  Thus  a  critical  feature  of  any  second- 
order-closure  model  is  how  this  macroscopic  scale  information 
is  folded  into  the  model.  Two  basic  approaches  are  currently 
used.  Either  the  macroscale  is  prescribed  over  the  flow  domain 
by  bulk  constraints,  or  a  modeled  differential  equation  is  used 
with  all  of  the  nonlocal  information  supplied  through  the 
prescribed  boundary  conditions.  In  either  case,  all  macroscale 
lengths  are  assumed  to  be  related  to  each  other  by  a  fixed  (for 
each  model)  ratio.  This  is  tantamount  to  assuming  that  the 
frequency  spectral  distributions  of  all  the  different  correla¬ 
tions  are  similar  when  properly  normalized.  This  assumption, 
of  course,  is  not  universally  valid. 
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The  shape  of  the  normalized  spectral  distributions  can 
have  significant  variations.  What  we  propose  herein  is  a 
modest  step  towards  allowing  the  model  to  use  shape  information. 
In  particular,  we  develop  a  model  in  the  next  section  which 
permits  the  ratio  of  the  length  scale  in  one  direction  to  the 
length  scale  in  the  direction  perpendicular  to  the  first  to 
vary  throughout  the  flow.  The  development  of  the  model  is 
followed  by  an  evaluation  of  the  new  model  coefficients  from 
data  obtained  in  the  atmospheric  unstable  mixed  layer.  This 
flow  is  chosen  because  of  the  strong  variation  in  the  ratio  of 
the  vertical  length  scale  to  the  horizontal  length  scale  which 
is  found  in  the  observations  of  the  surface  layer  in  such  a 
flow . 


> 
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2.  MODELING  CHOICES 


The  Reynolds  stress  equation  for  an  incompressible  fluid 
with  a  slightly  variable  density  in  the  presence  of  a  gravita¬ 
tional  field  may  be  written  rigorously  as 


Closure  of  the  system  of  equations  at  this  second-order  level 
requires  that  the  last  five  terms  in  this  equations  be  modeled. 
The  first  three  of  these  represent  turbulent  diffusion  terms. 

We  will  discuss  these  last.  We  believe  that  the  last  two  terms, 
the  pressure-strain  correlation  and  the  dissipation,  are 
generally  more  important. 

In  modeling  the  pressure-rate  of  strain  correlation, 

/  8u.  3u,\ 

=  p  f  >  we  aHow  for  tendency-towards-isotropy , 

as  proposed  by  Rotta  (1957)  and  generally  followed  in  second- 
order-closure  modeling,  on  the  broad  basis  that  pressure  fluc¬ 
tuations  randomize  directivity  since  the  pressure  is  a  scalar 
variable.  Several  refinements  to  Rotta's  tendency-towards- 
isotropy  model  have  been  proposed  in  the  literature.  These 
generally  involve  the  use  of  the  mean  strain.  We  have  not 
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included  such  a  term  in  our  considerations  here  for  two  reasons. 
First,  we  want  the  model  to  reduce  to  our  previously  developed 
single  scale  model  (Lewellen,  1977)  in  the  limit  of  an  isotropic 
scale  and,  second,  So  and  Peskin  (1980)  report  that  their  tests 
show  the  previously  suggested  refinements  to  be  of  questionable 
value  for  a  general  model. 

The  presence  of  a  wall  or  stratification  certainly  can 
introduce  a  bias  into  the  fluctuations.  This  directionality, 
which  can  be  characterized  for  our  purposes  by  the  normalized 
vector  (g^g^  =  1) ,  suggests  that  a  tendency  towards  axi- 

symmetry  about  can  also  be  included  in  modeling  the 

pressure-rate  of  strain  correlation.  The  most  general  symmetric, 
second-rank  tensor  which  can  be  formed  using  the  normal  or 
"vertical  vector"  is 


t 


ij 


+  BMj 


(2) 


where  A  and  B  are  arbitrary  scalars.  We  note  that  the  form 
given  for  t^j  can  also  be  deduced  from  the  assumption  that 
should  be  modeled  as  a  linear  combination  of  tendency- 
towards-isotropy  and  tendency-towards-axial  symmetry. 

The  tensor  that  we  wish  to  model,  P^  ,  has  zero  trace  as 
a  consequence  of  the  continuity  equation 


3u . 

^  -  0  <3> 

i 

and  thus  yields  no  net  contribution  to  the  energy  equation. 
Thus,  P^ j  represents  only  rearrangement  (or  scrambling)  among 
the  three  energy  components.  We  require  that  the  model  that  we 
choose  for  P^  should  have  the  same  property;  i.e.. 
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It  follows  that  we  have  the  relationship 

3A  +  B  =  0 


We  can  therefore  rewrite  t^  as 


cu  * 2A 


(i  sij  •  z  Mj) 


< 


(5) 


(6) 


A  factor  of  two  has  been  inserted  in  order  to  observe  the 
identity 


2A 


Si§j 
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The  tensor  is  the  projector  on  the  vertical  direction 

while  the  tensor  5^  -  g^gj  is  the  projector  on  the  horizon¬ 
tal  direction.  Our  tendency  model  thus  has  automatically  the 
property  that  a  source  or  sink  acts  on  the  vertical  component 
of  energy  depending  upon  the  sign  of  A  and  a  counterbalancing 
sink  or  source  acts  (isotropically)  on  the  two  horizontal 
components  of  energy. 

We  can  write  our  proposed  model  for  the  pressure-rate  of 
strain  correlation  as 
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To  complete  the  modeling  of  this  term,  we  must  relate  the  two 
time  scales  to  determinable  characteristics  of  the  turbulence. 
At  this  point,  we  also  introduce  the  constraint  that  our  two- 
scale  model  reduce  to  our  standard  single-scale  model  when  the 
horizontal  and  vertical  length  scales  are  equal.  At  the  other 
extreme  when  the  vertical  length  scale  is  much  smaller  than  the 
horizontal,  should  be  dominated  by  the  ratio  of  the 

vertical  length  scale  to  the  vertical  velocity  variance.  An 
expression  for  ,  which  provides  for  a  smooth  transition 
between  these  two  limiting  time  scales  is 


rs(l-Of+C.^  =  fFr 
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(10) 


This  expression  introduces  a  new  constant  C-^  .  The  exponent 
2/3  could  be  left  free  at  this  point  but  is  chosen  in  anticipa¬ 
tion  of  the  results  of  the  expansion  about  Ay/A^  -*-0  in  the 
next  section.  The  time  scale  for  the  tendency-towards-axial 
symmetry  should  depend  upon  how  anisotropic  the  components  of 


the  Reynolds  stress  are.  The  form  we  choose  is 
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which  introduces  two  more  constants,  5*  and  D  . 


Since  the  principal  contributions  to  dissipation  occur  as 
a  result  of  the  small  eddy  motions,  most  researchers  agree  that 
dissipation  should  be  modeled  as  an  isotropic  loss  term.  We 
add  to  this  the  requirement  of  a  smooth  transition  between 
appropriate  choices  for  the  two  limiting  cases  to  take  as  our 
choice 

v  /  5ui  9ui 
\3xk  3xk 

(12) 
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This  introduces  an  additional  constant,  C2  .  We  choose  to 
hold  b  fixed  at  1/8,  its  previous  standard  value.  A  simpler 
choice  for  dissipation  would  have  been  to  model  it  as  propor¬ 
tional  to  x^  ,  but  this  would  not  appear  to  give  the  correct 

dependence  on  a  in  the  limit  of  A  /A„  -*■  0  . 

w  v  H 

The  diffusion  terms  in  Eq.  (1)  are  combined  and  modeled 
simply  as 


<ukuiu j  > 
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^  ^  3u.u. 
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(13) 


We  have  refrained  from  introducing  any  new  model  coeffi¬ 
cients  in  the  diffusion  terms  and  have  stuck  to  the  philosophy 
of  a  simple  transition  from  the  previous  standard  model  when 
Av  *  A^  to  a  vertical  diffusion  which  depends  only  on  cw  and 
Av  when  Av  <<  A^  .  Models  which  attempt  to  incorporate  more 
of  the  physics  of  these  diffusion  terms  have  been  put  forward 
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by  others,  but  difficulty  with  the  accurate  determination  of 
the  required  model  coefficients  riste  any  real  gains  of  intro¬ 
ducing  additional  complexity  at  this  point. 

For  the  thermal  equations,  we  adopt  a  parallel  strategy. 
The  two  relevant  modeled  terms  are  the  "pressure  scrambling" 
of  the  heat  flux 


and  the  dissipation  of  the  temperature  variance 


(14) 


-  2k 


The  diffusion  terms  are  taken  as  in  Eq.  (13). 
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3.  EVALUATION  OF  NEW  MODEL  COEFFICIENTS 


The  new  formulation  given  in  the  last  section  introduces 
four  new  model  coefficients  C^,  C2,  £*  and  D.  These  coefficients 
will  be  determined  using  the  data  available  from  measurements 
in  the  unstable  surface  layer.  Spectral  distributions  taken  in 
the  surface  layer  under  such  unstable  conditions  clearly  show 
the  ratio  between  the  horizontal  and  the  vertical  scale  changing 
with  height.  Figure  1  from  Kaimal  (1973)  shows  tne  vertical 
scale  increasing  with  height  while  the  horizontal  scale  remains 
essentially  constant.  Data  taken  under  these  conditions  should 
provide  valid  information  for  determining  the  coefficients  in 
our  new  model. 


The  constant  in  the  dissipation  time  scale  can  be  evaluated 
by  considering  the  limit  of  free  convection.  In  this  limit  we 
expect  the  turbulent  kinetic  energy  equation  to  show  an  essential 
balance  between  buoyant  production  and  viscous  dissipation.  This 
assumes  that  the  remaining  diffusion  terms  are  less  important 
than  these  two  terms.  Data  as  analyzed  by  Wyngaard  (1979)  tend 


to  confirm  this  assumption.  This  balance  may  be  represented  as 

^  <1  -  e>  +  c2s  ^  at 


The  free  convection  limit  corresponds  to  z/L  -*■  -®.  Now  if 
we  also  consider  z/z^  ■*  0  then  we  can  expect  Av/A^  *  0  . 
Then  Eq.  (16)  reduces  to 

28^6%  _  2bC2  aw3 

To  =  A 


(17) 
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For  this  to  hold  when  A^/A  >>  1  (£  -*•  1) 

r  m  1.25  A  ,  0.65(1.25)  .  -  n 

C1  0779  z  ~-~~U7S9  2'° 


(24) 


It  thus  appears  appropriate  to  take  *  C2  =  2. 

The  remaining  coefficient  D  can  be  estimated  by  again 
appealing  to  the  free-convection,  surface-layer  limit  and 
considering  the  equation  for  any  one  of  the  normal  components 
of  the  Reynolds  stress.  To  the  same  approximation  as  Eq .  (16), 

the  equation  for  w'w*  may  be  written  as 
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When  this  is  combined  with  the  energy  equation  Eq .  (16), 
the  buoyancy  term  can  be  eliminated  to  yield 
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In  this  free  convection  surface  limit  we  have  already  seen 
that  £ 

be  further  reduced  to 


1,  -*■  A/2ow  and  ■+(A/2cP?)q^  ,  so  that  Eq .  (26)  can 
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Therefore  if  we  wish  to  permit  q  / 0^  to  become  large, it  is 

necessary  to  have  £*  =  1.  In  this  limit  it  is  necessary  to 

2  2 

keep  the  product  of  (1  -  £)q  /ow  in  evaluating  D. 

0  =  -o2(l  -  D)  +  q 2/ 3  ( 1  -  O  +  X  0 w 


(28) 
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If  we  take  *  0.2z^  then  the  two  separate  correlations 
of  Wvn/aard  for  cu/u  and  ;  /u  as vm:< 1 o t  i ca  1 1  v  lead  to 
(q- /c^)  (A/A  )~/3  =  3/4.  Thus  when  Eq .  (10)  is  used  for 
this  leads  to 


1  1 
'  5  "  e  =  rc 


Note  this  expansion  would  not  have  proceeded  snoot hlv  if 
some  other  exponent  rather  than  2/3  had  oeen  used  in  Eq.  (!'-•) 
No  new  coefficients  have  been  introduced  into  the 
ai f fusion  terms . 


ALGOR I  I'H'IS  FOR  SCALE  DE  i' F. A  v  I  NAT  I  ON 


The  use  of  separate  scales  for  the  horizontal  fluctua¬ 
tions  and  the  vertical  fluctuations  compounds  the  problem  of 


specifying  the  model.  The  single  isotropic  scale  use 


A  A  „ 


our 


previous  model  may  best  be  associated  with  the  vertical  scale 
in  the  present  model. 

Ratner  than  attempt  to  modify  the  scale  equation,  the 
equation  we  have  used  to  determine  /.  (Lewel1en,  1977)  will 
now  be  used  for  A... 
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Alternatively,  we  can  allow  the  vertical  scale  to  be  determined 
by  whichever  of  the  three  following  bounds  govern  locally: 


i) 


0.65 


ii )  <_  0 . 2  jz 


i  i  i  i 
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(32) 


The  first  of  these  recognises  that  changes  in  the  vertical 
sca.e  cannot  occur  abruptly.  When  this  first  condition  is 
combined  with  the  boundary  condition  that  the  vertical  scale  is 
Cere  i-rred lately  adjacent  to  a  surface,  it  yields  A  =  0.65  z  . 
The  second  condition  recognizes  that  the  vertical  scale  cannot 
exceed  a  certain  fraction  of  the  vertical  spread  of  the  region 
of  t  ..rbu enc  e  The  last  condition  recognizes  that  the  potential 
er.t  r  c  ••  ir.w  ived  with  the  vertical  fluctuations  in  a  stably 
stratified  fluid  cannot  exceed  the  kinetic  energy  in  the  fluc- 
tu.it  i,n  These  conditions  involve  only  a  slight  variation  from 
*  h.  st  used  f.  r  our  single  scale  model  (Lewellen,  1977). 

The  h,  r  i out  a  1  scale  cannot  be  specified  in  quite  so 
s  -  :  r  .  l.  •  rw  3rd  a  manner  It  must  not  satisfy  either  the  first 
ir  las:  c  f  the  conditions  in  Eq  .  (.32).  It  will  satisfy  the 

se  r.  d  -  nil  t  ion  if  the  two  scales  are  equal  in  the  center  of 
t :  t  f  1  -w  d  ,i  m  Sufficiently  close  to  the  surface,  the 
h  r  i  .*  -it  i  1  s .  ale  should  also  approach  zero  since  the  horizontal 
ve 1  >c it  v  flu.  ‘nation  is  constrained  to  approach  zero.  The  rate 
a'  w:  i  hr.  this  occurs  will  depend  on  the  surface  roughness  as 
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the  surface  shear  stress  serves  to  transfer  energy  from  the 
larger  scale  horizontal  fluctuation  to  the  small  scale  fluctua¬ 
tions  produced  by  the  large  vertical  gradients  in  the  larger 
scale  fluctuations.  Although  Figure  1  shows  remaining 

essentially  constant  while  A  decreases  with  decreasing  z  , 
there  is  a  gradual  increase  in  the  horizontal  energy  level  at 
high  frequency  corresponding  to  a  slight  decrease  in  average 
A^  with  decreasing  z  . 

Even  in  the  neutral  surface  layer,  there  can  be  consider¬ 
able  variation  in  how  A  /Au  behaves  close  to  a  surface.  Some 

V  n 

of  the  implications  of  this  can  be  seen  by  examining  the 
equilibrium  turbulent  kinetic  energy  equation  in  tne  neutral 
limit.  If  diffusion  is  neglected,  as  in  the  last  section,  this 
equation  may  be  written  in  terms  of  the  current  model  as 

2^L”  =  -  f(l  -  £)  +  2£  ^  3  1  (33) 

H  L  Av  q'J 

After  some  rearrangement,  this  may  be  rewritten  as 


(1  -  o 


(34) 


If  the  values  of  k  ,  b  ,  a  ,  and  °w/u*  used  in  evaluating 
Eqs .  (17)  -  (29)  are  used  here,  this  yields  an  equation  for 


(35) 


Thus,  for  A^/A^  varying  from  1  to  3 ,  this  leads  to  aw/q 
varying  from  0.53  down  to  0.3.  In  homogeneous  shear  flow  in  a 
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wind  tunnel,  Champagne  et  al .  (1970)  reported  a  value  of 
o  / q  =  0.49  while  Hinze  (1959)  reported  a  value  of  0.31  in 
turbulent  pipe  flow.  Values  between  these  have  been  reported 
for  other  wind  tunnel  and  atmospheric  studies.  These  differ¬ 
ences  are  in  the  right  direction  to  be  consistent  with  varia¬ 
tions  in  the  anisotropic  natures  of  the  eddy  scales.  That  is, 
A^/A^  should  be  closer  to  one  in  the  homogeneous  shear  flow 
case  than  for  pipe  flow. 

The  third  bound  of  Eq.  (32)  does  not  apply  to  A^  .  Rather 
in  a  stably  stratified  flow,  we  expect  the  occurrence  of  a 
natural  horizontal  scale  associated  with  the  ratio  of  the  mean 
horizontal  velocity  to  the  Brunt-Vai sala  frequency. 

We  are  not  proposing  a  specific  recipe  for  determining  A^ 
for  all  situations.  The  derivation  of  a  dynamic  equation  for 
A^  would  involve  considerations  similar  to  those  carried  out 
for  the  two-point  correlation  tensor  (Sandri,  1978).  This 
requires  considerably  more  work.  The  prescriptions  used  in  the 
next  section  will  demonstrate  the  type  of  specifications  we 
believe  are  required. 

5.  COMPARISON  WITH  EXPERIMENTAL  DATA 
We  will  first  compare  the  new  model  results  with  atmospheric 
surface  data.  This  will  be  done  using  a  superequilibrium 
version  of  the  model.  Then  complete  numerical  solutions  will  be 
used  to  compare  model  results  with  data  taken  across  the  full 
height  of  the  unstable  free  convection  layer. 
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a )  Surface  Layer  Comparisons 

When  the  turbulent  layer  immediately  adjacent  to  a  surface 
is  considered,  the  time  scale  of  the  turbulence  is  sufficiently 
small  that  the  ensemble  mean  values  can  be  assumed  to  be  at 
their  steady-state  values.  When  this  condition  is  combined 
with  the  assumption  of  horizontal  spatial  homogeneity  and  high 
Reynolds  number,  the  momentum  and  energy  equations  reduce  to 
the  conditions  that  the  vertical  momentum  and  heat  flux  remain 
constant  with  respect  to  variations  in  height.  The  Reynolds 
stress  equations  can  then  be  solved  for  the  relationships 
between  the  turbulent  correlations  and  the  mean  velocity  and 
temperature  gradients.  This  can  be  done  with  the  turbulent 
diffusion  terms  included  but,  since  our  model  plays  down  the 
relative  role  of  these  diffusion  terms,  we  choose  to  first 
compare  the  model  results  with  the  data  with  these  diffusion 
terms  neglected.  This  constitutes  what  Donaldson  (1973)  has 
termed  the  superequilibrium  approximation. 

When  the  superequilibrium  approximation  is  applied,  Eq. 

(1)  reduces  to 
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To  complete  the  system  of  equations,  we  must  also  add  the 
heat  flux  and  temperature  variance  equations: 
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(37) 


with  y  arbitrarily  chosen  to  provide  a  relatively  rapid 
transition  between  a  neutral  value  equal  to  az  for  very 
small  z  and  the  constant  value  proportional  to  the  mixed  layer 
height  at  larger  z  .  As  argued  in  the  last  section,  the  actual 
transition,  i.e.,  y  ,  should  depend  on  the  surface  roughness 
in  any  given  case.  On  the  stable  side  (z/L  >0),  av  is  equal 
to  az  or  0.2L,  whichever  is  smaller,  but  no  limit  is  placed 
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on  Au  .  It  is  permitted  to  grow  proportional  to  z  for  all 
z/L  >  0  .  Again,  the  more  precise  variation  of  A^  should 
depend  on  surface  roughness  with  a  transition  from  the  surface 
small  scale  eddies  to  the  larger  internal  wavelengths  associated 
with  the  ratio  of  the  mean  horizontal  velocity  to  the  Brunt- 
Vaisala  frequency.  Even  after  the  transition,  A^  will  continue 
to  grow  proportional  to  z  because,  as  will  be  seen,  9u/3z 
is  constant  with  respect  to  variations  in  z/L  after  the 
critical  Richardson  number  is  reached. 

Figures  2,  3,  and  4  show  that  there  is  not  much  variation 
between  the  previous  single-scale  model  results  for  $  ,  or 

ow/u^.  and  the  results  for  the  current  two-scale  model.  Either 
provides  a  creditable  representation  of  the  data.  The  principal 
difference  shows  up  in  the  horizontal  variances.  The  single¬ 
scale  model  predicts  a  variation  with  respect  to  z/L  only 
which  is  not  observed  in  the  data.  Rather  the  data  are  found 
to  be  relatively  independent  of  z/L  and  to  depend  on  the  mixed 
layer  height  instead.  Figure  5  shows  the  dependence  of 
oH/u*  ,  °h  =  ?  (°w  +  °v  )  ]  on  Zi'/L  anc3  Z/L  as  Predicted 
by  the  present  model.  The  single-scale  model  would  show  no 
dependence  on  z^/L  . 

In  the  present  model  modification,  we  have  not  adjusted 
any  coefficients  in  the  temperature  correlation  equations.  How¬ 
ever,  Figures  6  and  7  show  that  the  modification  does  provide 

a  modest  improvement  in  representing  the  data  for  ue/u*e*  and 
~2  2 

6  /e*  .  The  peak  in  these  correlations,  which  we  previously 


O 
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Figure  7.  Dimensionless  rms  temperature  fluctuations  as  a 

function  of  the  Monin-Obukhov  similarity  variable 

(Lewellen  and  Teske,  1973)  single  scale  model; 

-  two-scale  model;  data  from  Wyngaard,  et  al . 

(1»71)] 
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felt  was  due  to  experimental  error  near  z/L  =  0  where  6^=0, 
is  partially  reproduced  with  the  present  model,  although  it  is 
shifted  towards  positive  z/L  . 

b )  Unstable  Mixe d  Layer  Comparison 

The  next  test  of  the  two-scale  model  was  to  repeat  the 
calculations  of  turbulence  in  an  unstable  mixed  layer  (Lewellen, 
Teske  and  Donaldson,  1976)  using  a  numerical  code  which  computes 
the  full  second-order  closure  set  as  a  function  of  z  and  t  . 
Horizontal  homogeneity  is  assumed.  For  this  free  convection 
problem,  a  constant  positive  heat  flux  is  applied  at  the  surface 
and  a  stable  temperature  lapse  rate  is  assumed  at  the  top  of  the 
domain.  With  the  mean  velocity  set  to  zero,  a  free  convection 
mixed  layer  forms  adjacent  to  the  surface.  The  thickness  of 
this  mixed  layer  increases  with  time,  but  the  velocity  variances 
exhibit  similarity  as  a  function  of  z/z^  when  normalized  by 
the  characteristic  velocity 

«*  -  ^  zi>/T0]1/3  <41) 

with  z^  the  depth  of  the  mixed  layer.  Figures  8,  9  and  10 
compare  the  model  results  with  data  from  Willis  and  Deardorff 
(1974) .  Both  models  represent  the  vertical  velocity  variance 
very  well.  The  major  difference  is  in  the  horizontal  velocity 
variance.  The  two-scale  model  provides  a  much  better  represen¬ 
tation  of  near  the  surface  and  some  improvement  near  the 

inversion.  It  also  provides  improvement  in  the  representation 
of  the  temperature  variance  near  the  inversion. 
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Figure  8.  Comparison  of  free  convection  prediction  and  experi¬ 
ment  for  the  normalized  vertical  velocity  fluctuations 
WW  =  wTvTr/w*:  ,  as  a  function  of  normalized  height 
Z  =  z/z-y  ,  where  w*  is  the  characteristic  convective 
velocity  and  z^_  is  the  depth  of  the  mixed  layer. 

Model  predictions  are  for  the  two-scale  model  ( - ) 

and  the  single-scale  model  ( - ).  The  solid  data 

points  are  for  laboratory  observations  given  by 
Willis  and  Deardorff  (1974). 


6 .  CONCLUDING  REMARKS 


Although  the  comparisons  between  model  and  data  shown  here 
suggest  some  further  adjustments  in  scale  modeling  and/or  more 
sophisticated  diffusion  modeling  are  desirable,  they  do  indicate 
that  the  use  of  a  two-scale  model  has  the  potential  for  signifi¬ 
cantly  improving  model  predictions.  This  is  particularly  true 
for  the  correlations  involving  fluctuations  parallel  to  a 
surface . 

Other  applications  of  turbulence  closure  models  to  the  free 
convection  problem  may  be  found  in  the  work  of  Mellcr  and 
Durkin  (1975),  Garwood  (1977),  Zeman  and  Lumley  (1976),  Lumley, 
Zeman  and  Siess  (1978)  ,  Wyngaard  (1979) ,  and  Warn-Varmas  and 
Piacsek  (1979)  .  The  last  four  of  these  attempt  to  improve  upon 
the  modeling  by  bringing  in  more  information  about  the  triple 
velocity  correlation.  This  assumes  that  the  basic  difficulties 
in  the  modeling  are  due  to  inadequate  representation  of  the 
local  single-point  higher  correlations  on  the  local  value  of 
the  lower  order  correlations  and  their  derivatives.  We  have 
chosen  the  alternati/e  approach  of  bringing  in  more  information 
about  the  macroscale  flow.  Due  to  the  existence  of  large  eddies 
in  the  flow,  we  believe  it  is  more  important  to  bring  in  some 
rough  information  about  the  nature  of  the  two-point  correlations 
than  it  is  to  delve  more  deeply  into  the  intricacies  of  the 
relationships  between  the  single-point  higher  order  correla¬ 
tions  and  the  derivatives  of  the  single-point  lower  order 
correlat ions . 
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ABSTRACT 

The  interdependence  between  the  transport  of  water  vapor 
through  the  atmospheric  boundary  layer  and  its  distribution 
within  the  layer  is  reviewed.  The  problem  is  approached  using 
second-order  closure,  turbulent  transport  modeling  with  the 
results  related  to  simple  mass  transfer  coefficients.  The 
dependence  of  moisture  fluctuations  on  such  parameters  as 
wind  speed,  thermal  stability,  and  surface  roughness  is 
summarized.  The  results  of  example  calculations  are  used  to 
illustrate  such  phenomena  as:  the  influence  of  water  vapor 
flux  on  boundary  layer  stability;  the  interaction  of  water 
vapor  with  a  strong  temperature  inversion  at  the  top  of  the 
boundary  layer;  the  interaction  between  thermal  radiation, 
condensation,  and  turbulence  when  phase  change  occurs  in  the 
boundary  layer;  and  the  relative  role  of  different  size  eddies 
in  transporting  water  vapor  across  the  boundary  layer. 


This  work  was  supported  in  part  by  the  Naval  Air  Systems  Command, 
Contract  No.  N00019-79-C-0366 . 
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INTRODUCTION 


The  atmospheric  boundary  layer,  which  consists  of  approxi¬ 
mately  the  lowest  1  km  of  the  atmosphere,  plays  a  key  role  in 
controlling  the  cycling  of  water  thru  the  atmosphere.  On  a 
global  average  basis,  the  turbulent  transport  of  water  vapor 
thru  this  atmospheric  layer  must  balance  the  average  precipi¬ 
tation  of  liquid  water,  about  1  m  per  year  (Coantic,  1978). 

This  transport  is  achieved  by  turbulent  eddies  mixing  the  air 
with  high  relative  humidity  next  to  the  surface  with  dryer  air 
above  it.  Thus  the  distribution  of  water  vapor  within  the 
boundary  layer  is  inherently  dependent  upon  the  transport  of 
water  vapor  thru  it.  It  is  equally  dependent  upon  the  dynamics 
of  the  turbulent  air  motion  in  this  layer. 

Although  it  is  the  turbulent  eddies  which  drive  the  trans¬ 
port  of  vapor,  it  is  just  as  valid  to  note  that  the  latent 
energy  transported  thru  the  boundary  layer  with  the  water  vapor 
is  ultimately  responsible  for  fueling  most  of  the  atmospheric 
motions  producing  the  turbulence.  As  pointed  out  in  the  previous 
paper  by  Klemp,  the  intensity  of  an  atmospheric  storm  can  depend 
critically  upon  the  water  vapor  content  of  the  boundary  layer  air 
which  feeds  into  the  storm. 

To  understand  the  transport  of  water  vapor  it  is  necessary 
to  understand  the  dynamics  of  the  turbulent  eddies  ranging  in 
scale  from  approximately  10  ym  to  100  m  next  to  the  surface. 

The  smallest  eddies  dominate  the  transport  close  to  the  surface. 
At  increasing  altitudes,  the  scale  of  the  eddies  most  responsible 
for  vapor  transport  will  also  increase.  Since  we  are  limiting 
our  attention  to  a  boundary  layer  thickness  of  order  1  km  or 
less,  the  largest  scale  eddy  of  interest  is  of  order  of  a  few 
100  m.  Larger  scale  motions  must  be  broken  down  to  smaller 
scales  by  shear  production  before  it  enters  directly  into  the 
vertical  transport  of  vapor  close  to  the  surface. 
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In  the  next  section,  a  review  of  turbulent  transport  modeling 
is  given  using  a  second-order  closure  approach.  This  allows  some 
of  the  physics  of  the  production  and  decay  of  turbulence  to  be 
incorporated  into  a  semi-empirical  framework  for  investigating  the 
sensitivity  of  vapor  transport  to  such  variables  as  wind  speed, 
thermal  stability,  surface  roughness,  etc.  These  turbulent  trans¬ 
port  equations  can  also  be  integrated  to  indicate  to  which  para¬ 
meters  the  simple  bulk  mass  transfer  coefficient  is  sensitive. 
Results  of  several  example  calculations  are  reviewed  in  later 
sections . 
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MODEL  EQUATIONS 


The  ensemble-averaged,  Eulerian  equation  of  mass  continuity 
for  a  species  concentration  may  be  written  as 


DC 

Dt 


+  D 


(1) 


where  C  is  the  mean  value  of  the  concentration  of  a  species  such  as 
water  vapor,  c  is  the  fluctuating  value,  u.  is  the  fluctuating 
velocity  in  direction  ,  D  is  the  molecular  diffusion  coefficient 
and  S  any  source  term.  Equation  (1)  is  exact  but  undetermined  even 
if  the  mean  velocity  is  known  because  of  the  presence  of  the  addi¬ 
tional  variable  uTc  representing  the  ensemble  average  of  the  corre¬ 
lation  between  the  fluctuating  velocity  u^  and  concentration 
fluctuations  at  a  point.  The  vertical  component  of  this  turbu¬ 
lent  transport  of  species  is  the  primary  variable  of  interest  to 
us  here. 


An  exact  equation  can  be  written  for  this  variable  (e.g., 
Donaldson ,  1973) . 
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The  first  two  terms  on  the  right-hand  side  are  production  terms 
due  to  the  interaction  of  correlations  and  the  mean  concentration 
or  velocity  gradients.  The  third  term  is  a  buoyancy  term  repre¬ 
senting  the  direct  effect  of  a  correlation  between  the  species 
and  virtual  potential  temperature  fluctuation  0  .  The  co¬ 

efficient  to  this  latter  term  is  the  ratio  of  the  gravitational 
acceleration  g  to  the  reference  temperature  Tq  . 

The  last  four  terms  on  the  right-hand  side  introduce  variables 
other  than  second-order  correlations  and  they  leave  the  system  of 
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equations  undetermined.  The  task  of  second-order  closure  is  to 
model  these  terms  as  functions  of  the  second-order  correlations 
and  the  mean  flow  variables.  A  simple  modeled  form  appropriate 
for  high  Reynolds  number  flow  and  providing,  at  least,  the  minimum 
information  needed  to  close  the  system  may  be  written  as 
(Lewellen ,  1977) 
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The  coefficients  to  the  last  two  terms  are  carried  over  from 
modeling  the  velocity  and  temperature  correlation  equations,  while 
q  and  A  represent  velocity  and  length  scales  of  the  turbulence. 
When  dealing  with  a  plume,  it  appears  necessary  to  recognize  that 
the  velocity  and  length  scales  for  the  species  fluctuations  can 
be  quite  different  from  those  of  the  velocity  or  temperature 
fluctuations  (Lewellen  and  Teske,  1976),  but  for  the  spatially 
homogeneous  case  of  most  interest  here,  this  problem  will  not 
arise . 

We  do  not  expect  that  the  last  two  modeled  terms  in  Eq .  (3) 
used  to  replace  the  complex  terms  of  the  exact  equation  will 
faithfully  represent  all  of  the  information  present.  However, 
for  most  problems,  we  are  interested  in  only  a  small  part  of  the 
information  contained  in  the  complete  turbulent  spectrum.  We  be¬ 
lieve  that  the  two  modeled  terms  provide  at  least  the  minimum 
amount  of  desired  information  needed  to  close  the  system  at  the 
second  order.  The  first  modeled  term  introduces  diffusion  to 
prevent  excessive  gradients  in  the  species  flux.  The  other 
modeled  term,  a  tendency- towards-isotropy  term,  introduces  the 
required  feedback  which  permits  the  flux  to  reach  an  equilibrium 
level  even  in  the  presence  of  large  production  contributed  by  the 
first  three  exact  terms  on  the  right-hand  side  of  Eq .  (3). 

The  effect  of  buoyancy  on  transport  comes  into  Eq .  (3)  in 
two  ways:  through  its  influence  on  the  stability  of  the  velocity 
fluctuations,  and  through  the  buoyant  term  involving  the  cross 
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correlation  of  the  species  and  the  virtual  potential  temperature,  c6 
which  appears  directly  in  Eq.  (3).  This  term  is  not  a  result  of 
our  closure  modeling  but  arises  directly  from  the  buoyant  term 
in  the  momentum  equation.  However,  modeled  terms  must  appear  in 
the  equation  derived  for  c6  .  If  these  are  treated  in  a 
similar  fashion  to  those  in  Eq.  (3),  the  equation  for  c6v  may 
be  written  as 


Dcev 

Dt 


“  3C_ 

v  ax. 
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+  0.3 


ace 


O.A5qc6v 
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With  the  velocity  and  temperature  fields  specified  or  calcu¬ 
lated  from  similarly  modeled  equations  (Lewellen,  1977),  Eqs .  (1)  , 
(3) ,  and  (4)  form  a  closed  set  to  determine  C  ,  wc  ,  and  c 0 
when  appropriate  boundary  conditions  are  applied. 

The  humidity  transport  can  influence  the  velocity  and 
temperature  fluctuations  even  in  the  absence  of  phase  change 
through  the  buoyancy  terms  involving  0v  since 


9  =  T  -  T  +  ^+  0.61TH=(? 

V  o  c  o 
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with  H  the  mean  water  vapor  mixing  ratio,  c^  the  specific 
heat  of  air,  and  P  the  mean  potential  temperature.  In  the 
marine  atmospheric  boundary  layer,  the  contribution  of  the 
moisture  flux  to  buoyancy  can  sometimes  exceed  the  contribu¬ 
tion  from  the  temperature  difference.  In  terms  of  the  Bowen 
ratio  (£  =  CpW?/Lwh  with  JZ  the  latent  heat  of  evaporation) 

—5 —  —5-/1  1  0.61  c  p^  o  \ 
wev  -  we  (l  +  — g —  — ) 


V0 


we  (l  + 


Typical  values  of  6  for  the  air-sea  interaction  are  =0.05. 
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INFLUENCE  OF  STABILITY  IN  THE  SURFACE  LAYER 


For  the  steady  state  spatially  homogeneous  case  which  should 
be  appropriate  close  to  the  surface,  but  above  any  sources  or 
sinks  of  the  surface  elements  themselves,  Eq .  (1),  (3),  and  (4) 
reduce  to 


wc  =  const. 
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When  the  diffusion  term  in  Eq .  (7)  is  neglected  on  the  basis  that 
c {  is  approximately  constant,  then  Eqs .  (5) -(7)  can  be  mani¬ 
pulated  to  give  r  —  -i 

tiS  +  ^  __Auev 

5S  ,  .  ?.3.  »533_  Aq  (8) 

1  +  S_  _12_ 

L  6o  0.34?  ?2J 

The  term  in  the  bracket  is  a  function  of  the  Monin-Obukhov 

3  _  2  _  _ 

similarity  variable  z/L  (L  =  -T0u*/kgwcv  with  u* =  uwq + vwq) 
only.  In  terms  of  the  usual  similarity  variable  bc  ,  Eq .  8  can 
after  some  manipulation  be  written  as 

u,.kz  3C  +  ,  ...  -3kzu*q/4w 

-r_-  (z/L)  = - - - 3 - - —  (9) 

wc  02  c  [1  -  11.3(u;/qww) (z/L) ] 

with  k  ,  von  Karman's  constant,  equal  to  the  value  necessary  to 
make 

q  (0)  =  —  ~  =1 

^mv  '  Uj.  3z 


In  Eq .  (9)  3fv/3z  has  been  replaced  by  (we^/wc) jC/ ^z  in 
recognition  of  the  fact  the  model  treats  <t>c  as  identical  to 
i ~  .  This  function  as  obtained  from  the  temperature  data  of 

Businger,  et.  al.,  1971  is  shown  in  Figure  1,  along  with  the 
predictions  using  a  second-order  closure  model  with  the  cif- 
fusion  terms  included  (Lewellen  and  Teske,  1973).  In  this 
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calculation,  A  was  taken  equal  to  0.65z  except  on  the  stable 
side  where  it  has  an  upper  bound  equal  to  0.2L  .  Most  data 
appear  to  agree  that  =  'f'e  >  although  some  data  have  been  re¬ 
ported  at  variance  with  this  theoretical  result  (e.g.,  Verma, 
Rosenberg  and  Bladd,  1978).  Two  factors  which  can  sometimes  allow 
che  humidity  to  be  transported  by  the  eddies  in  a  different  manner 
from  the  temperature  are  thermal  radiation  at  sufficiently  low 
wind  speeds  and  the  possibility  of  the  boundary  conditions  being 
introduced  at  different  scales.  These  factors  should  be  the  ex¬ 
ception  rather  than  the  rule. 


Figure  1  shows  that  humidity  gradients  will  be  strongest 
where  there  is  either  strong  transport  of  humidity,  low  winds, 
or  where  the  flow  is  quite  stable. 


Humidity  transport  data  are  often  reported  in  the  form  of 
a  bulk  transfer  coefficient 


C 

q 


(10) 


with  l’r  and  Hr  taken  at  some  arbitrary  reference  height  such 
as  at  10m  .  Equation  (9),  in  combination  with  a  similar  but 
quantitatively  different  function  for  SU/3z  ,  can  be  inte¬ 

grated  from  the  effective  surface  roughness  height  zq  to  z r  to 
yield  as  a  function  zc/zr  an^  zr/L  .  However,  attempts 

at  correlating  data  in  this  fashion  have  only  met  with  modest 
success  as  illustrated  in  Figure  2  fro  Smith,  1974.  The  rela¬ 
tively  large  scatter  in  the  data  is  not  adequately  explained  by 

the  uncertainty  involved  in  the  effective  z  for  a  water  sur- 

o 

face  with  waves.  It  probably  also  represents  some  difficulty  in 
measuring  the  true  surface  temperature  of  the  water. 

The  problem  of  determining  the  bulk  transfer  coefficient 
for  transport  across  the  total  boundary  layer  can  be  subdivided 
into  the  problem  of  transporting  across  the  separate  layers 
which  make  up  the  boundary  layer.  There  are  at  least  four 
separate  regions  which  can  usefully  be  defined:  1)  the  surface 
sublayer  or  canopy  layer,  2)  the  surface  constant  flux  layer, 
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3)  the  main  boundary 
layer.  Equation  (10) 
layer  as 


layer,  and  4)  any 
can  be  rewritten 

C  =  (U  I.R.)'1 
q  r  i  l 


capping 
for  the 


stable  inversion 
total  boundary 

(ID 


where  R.  is  defined  as  an  individual  laver  resistance 
l 

R.  =  AH./wh"  (12) 

1  1  s 

The  best  defined  region  is  Region  2  where  R9  is  a 
function  of  thickness  and  L  only.  For  a  neutral  layer  it  is 
given  simply  by 

Uj_R2  =  1.9  £n(z./Zg  (13) 


where  Zg  ^  is  the  top  of  Region  1.  The  influence  of  stability 
on  R2  through  the  Monin-Obukhov  length  L  is  given  in  Figure  3 

The  sublayer  "Region  1"  resistance  is  usually  parameterized 
in  terms  of  an  effective  roughness  zq  .  The  effective  roughness 
for  momentum  will  not  be  equal  to  the  effective  roughness  for 
mass  transfer,  in  general. 


To  better  understand  the  difference  between  the  effective 
roughneos  for  momentum  and  the  effective  roughness  for  moisture 
transport,  it  is  necessary  to  look  at  the  transition  layer  which 
generally  consists  of  some  canopy  of  vegetation  over  land  and  a 
wavy  surface  over  water.  Integration  of  the  momentum  and  species 
equation  across  this  transition  or  canopy  layer  of  height  h  with 
appropriate  source  and  sink  terms  proportional  to  the  wetted  area 
per  unit  volume,  ,  leads  to 


(14) 
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where  D/ v  is  the  ratio  of  the  molecular  diffusivity  of  water 
vapor  in  air  to  the  kinematic  viscosity  of  air,  D  /D^  is  the 
ratio  of  pressure-to-friction  forces  on  the  surface  elements. 

The  1/4-power  exponent  appearing  in  Eq.  (14)  is  only  appropriate 
for  D/ v  >  1  .  The  ratio  of  the  two  integrals  is  a  function  of 
element  shape,  S  and  Reynolds  number,  Re  so  that 
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Under  standard  conditions,  D/v  r  1.5  for  water  vapor  diffusing 
in  air,  so  for  a  rough  approximation  setting  G  =  1  , 
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(16) 


The  ratio  of  pressure  drag  to  friction  drag  should  be  ex¬ 
pected  to  increase  with  wind  until  wave  breaking  occurs  to 
greatly  increase  the  contact  surface  area  between  air  and  sea. 

Thus,  C  /CL  should  first  decrease  with  wind  speed  and  then 

q  D  r 

increase  as  wave  breaking  occurs. 

Within  a  vegetation  canopy  the  drag  ratio  depends  upon  the 
structure  of  the  canopy.  It  has  values  of  3  or  greater  for  a 
typical  forest  canopy  (Thom,  1975) . 

Examples  of  calculations  of  the  resistance  across  the  outer 
boundary  layer  and  the  upper  inversion  layer  will  be  given  in 
the  following,  sections. 

In  considering  transport  thru  the  main  boundary  layer  we 
will  concentrate  on  conditions  leading  to  an  unstable  mixed  layer 
capped  by  a  stable  temperature  inversion.  Three  reasons  for 
this  biased  treatment  are:  1)  that  on  the  average  the  boundary 
layer  is  unstable;  2)  the  majority  of  water  vapor  transport  occurs 
under  unstable  conditions  since  R2  is  sufficiently  large  under 
stable  conditions  to  impede  transport;  and  3)  it  is  simpler  than 
for  stable  conditions,  The  reader  interested  in  second-order 
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closure  model  results  of  the  inherently  unsteady,  stable  layer 
is  referred  to  Lewellen  et  al.  1974,  Yamada  and  Mellor  1975, 
or  Borst  and  Wyngaard  1978.  Simple  parameterized  models  of 
the  stable  layer  have  been  recently  given  by  Zeman,  1979  and 
Yamada,  1979. 
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TRANSPORT  THRU  THE  UNSTABLE  MIXED  LAYER 


Typical  profiles  of  potential  temperature  and  water 
mixing  ratio  in  an  unstable  mixed  layer  are  given  in  Fig.  4. 

The  data  shown  here  were  taken  from  the  Air  Mass  Transportation 
Experiment  (AMTEX)  as  reported  by  Wyngaard,  et  al.  1978.  The 
surface  layer  relations  hold  to  altitudes  of  approximately 
0.2  and  then  the  layer  is  essentially  uniformly  mixed  up 
thru  approximately  0.9  Z.  .  The  moisture  content  in  the  mixed 
layer  is  a  function  of  time  but  remains  a  function  of  z/L 
and  zq/L  only. 

The  top  portion  of  the  mixed  layer  and  the  capping 
inversion  layer  introduce  additional  parameters  into  the 
problem,  namely  temperature,  humidity,  and  velocity  jumps 
across  this  thin  sharp  gradient  region.  As  long  as  phase 
change  does  not  occur,  the  similarity  between  temperature  and 
humidity  transport  should  remain  valid  in  the  form 


Kw_  _  AH 

0w  AO 


(17) 


Since  AH/A9  is  generally  negative  the  ratio  of  the  fluxes 
must  also  be  negative.  Figures  5  and  6  show  data  on  the 
moisture  and  heat  flux  variation  within  the  mixed  layer. 

The  complementary  temperature-humidity  correlation  is  given 
in  Fig.  7.  The  temperature-humidity  correlation  is  positive 
at  the  bottom  and  negative  at  the  top.  The  curves  represent 
the  results  of  A. R. A. P.'s  model  simulation  of  correlations 
approximating  those  reported. 

The  available  data  does  not  meet  the  conditions  necessary 
for  a  rigorous  model  verification  test  (Lewellen  et  al.  1979). 
The  model  is  started  with  the  mean  flow  variables  approximating 
those  reported  for  15  February  and  the  turbulent  fluctuations 
allowed  to  adjust  to  these  mean  distributions.  The  geostrophic 
wind  was  adjusted  to  give  the  reported  value  of  the  surface 


Figure  6.  Vortical  humidity  flux  as  a  function  of  attitude  for  conditions 
simulating  the  Air  Mass  Transformation  Experiment  (AMTKX) . 
o-data,  Wyngaard  et  al.  1978. 


shear  stress  and  the  sea  surface  temperature  was  adjusted 
to  yield  the  reported  value  of  surface  heat  flux.  These 
adjustments  also  fix  the  surface  humidity  flux  which  agrees 
with  the  reported  value  to  within  10  percent. 

In  the  lower  half  of  the  mixed  layer,  the  model  predictions 
fall  well  within  the  scatter  of  the  data.  There  are  marked 
disparities  in  the  upper  half  of  the  mixed  layer,  but  there 
appears  to  be  a  good  reason  for  this  disparity.  The  model 
predicts  relatively  dense  scattered  clouds  between  800  m  and 
1355  m  (maximum  liquid  water  content  equal  to  0.69  gm/kg)  and 
the  authors  of  Ref.  23  report  stratocumulus  clouds  over  as 
much  as  80  percent  of  the  area  in  the  upper  five  hundred 
meters  of  the  mixed  layer.  Since  cloud  water  affected  their 
temperature  and  humidity  measurements,  they  attempted  to 
stay  out  of  the  clouds  while  taking  the  data.  Thus  the  data 
points  represent  the  profiles  in  the  cloud-free  environment 
between  the  clouds  while  the  prediction  line  represents  an 
ensemble  average  including  the  clouds. 

When  there  are  no  clouds  present  hw  is  positive  at 
while  wc  is  negative  as  predicted  by  Eq .  17.  However 
when  clouds  form  in  the  upper  portions  of  the  mixed  layer 
thermal  radiation  becomes  important  and  destroys  the  similarity 
between  temperature  and  humidity  profiles.  The  curves  shown  in 
Figs.  4-6  include  water  vapor  condensation  and  thermal  radiation 
as  detailed  by  Oliver  et  al.  1978.  The  short  absorption  depth 
of  the  infrared  radiation  within  a  cloud  ensures  that  cooling 
is  confined  to  a  relatively  thin  layer  next  to  the  top  of 
the  cloud.  Turbulent  transport  delivers  the  heat  from  down  in 
the  cloud  to  the  cloud  top  to  balance  the  radiational  cooling 
at  the  top.  Thus  in  the  presence  of  a  cloud  Kw/ew  is  positive 
at  even  when  AC/ AH  <  0  . 


BOUNDARY  LAYER  WITH  STRATUS  CLOUD 


A  boundary  layer  over  a  large  expanse  of  water  capped 
by  a  subsidence  inversion  will  turbulentlv  ingest  moisture 
from  the  surface.  Long-waveleneth  cooling  of  the  water-vapor¬ 
laden  air  may  then  lead  to  condensation  somewhere  within  the 
layer.  In  an  expansive  basin  or  region  influenced  only  by 
adveetion  in  the  form  of  a  steady  subsidence,  a  periodic  steadv 
state  may  form  in  which  a  solar  driven  diurnal  cycle  exists 
about  a  mean  boundary  layer  state  with  a  steady  mean  inversion 
height  controlled  by  the  radiative  conditions  (season  and 
latitude),  surface  temperature,  and  upper  conditions  (wind, 
water  content  and  subsidence)  set  by  the  larger  circulation. 

Such  an  episode  was  illustrated  by  Oliver,  et  al .  1978  for 
conditions  approximating  the  summer  stratus  layer  along  and 
off  the  coast  of  California. 

The  conditions  for  this  illustration  were  selected  to 

generally  conform  to  the  conditions  described  by  Neiberger 

(19  w)  and  more  recently  observed  by  Mack  et  al  .  (1974). 

Surface  boundary  conditions  are  established  for  a  fixed 

temperature  sea  state  of  17°C,  surface  water  mixing  ratio 

at  saturation,  and  surface  roughness  determined  by  Froude 

scaling.  Radiation  conditions  are  for  latitude  40°N  at 

summer  solstice.  Geostrophic  conditions  set  an  upper  level 

wind  ,  U  ,  of  10  m  s  ^  and  subsidence  whose  characteristic 
£  -  5  - 1 

value  is  0.5  x  10  s  (Neiberger  et  al . ,  1961;  Lilly,  1968). 
Initial  conditions  (which  will  be  lost  after  several  days  of 
simulation)  were  selected  to  correspond  to  a  mildly  stable 
temperature  profile  of  lapse  rate  3?v/?z  =  0 . 003°C  m'1  and 
clear  skv  with  an  initial  relative  humidity  of  0.9.  The  cal¬ 
culation  is  begun  shortly  before  sunset  the  first  day  and  runs 
for  six  days.  Behavior  during  daylight  hours  must  be  considered 
somewhat  approximate  in  the  present  illustration  because 
droplet  scattering  of  direct  solar  radiation  would  reduce 
solar  penetration  into  the  cloud  interior. 
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The  quasi-periodic  evolution  of  cloudiness  in  the 
boundary  layer  is  shown  in  Fig.  8.  The  stratus  grows  both 
downward  and  upward  during  nocturnal  periods  and  thickens 
until  sunrise  when  the  solar  ht  ting  begins  to  disperse  the 
stratus.  Radiative  cooling  is  dominant  at  the  cloud  top  and 
produces  an  unstable  lapse  rate  within  the  cloud.  A  cor¬ 
respondingly  enhanced  turbulence  production  is  also  maximum 
in  nocturnal  periods  and  in  the  upper  portion  of  the  cloud. 

The  radiation  cooling  at  the  cloud  top  strengthens  the  capping 
inversion  which  oscillates  up  and  down  over  the  period  of  the 
diurnal  cycle. 

We  note  that  the  predicted  cloud  base  in  these  diurnal 
cycles  is  highest  in  late  afternoon  and  lowest  in  early 
morning  (Fig.  8)  while  cloud  top  is  highest  in  early  morning 
and  lowest  in  late  afternoon.  Correspondingly  the  inversion 
heicht  is  highest  in  early  morning  and  lowest  in  later  after¬ 
noon  -  a  result  which  is  opposite  to  that  for  inversion  height 
cycles  driven  by  solar  heating  in  the  cloud-free  boundary 
layer  (Lewellen  et  al . ,  1974).  This  predicted  result  is  in 
accord  with  the  stratus  observations  of  Neiberger  (1944)  as 
well  as  those  of  Mack  et  al,  (1974)  in  which  they  detected 
a  diurnal  variation  of  inversion  height  in  the  presence  of 
California  coastal  stratus  which  regularly  showed  maxima  in 
the  early  morning.  It  is,  of  course,  the  r ad i a t i ve- turbul en t 
drive  of  cloud-top  cooling  which  allows  the  radiating  layer  to 
propagate  condensation  upward  during  the  night  while  turbulence 
cools  the  cloud  interior  below. 


Cloud  bo»e,  h,  or  top 


RELATIVE  CONTRIBUTION  OF  DIFFERENT  SIZE  EDDIES 


The  unstable  boundary  layer  can  often  get  organized  into 
eddies  which  extend  the  full  height  of  the  mixed  layer  in 
the  form  of  longitudinal  roll  vortices  aligned  close  to  the 
direction  of  the  mean  wind.  An  example  of  the  clouds  resulting 
from  this  motion  are  shown  in  Fig.  9.  When  the  boundary  layer 
is  treated  as  spatially  homogeneous  then  the  motion  of  these 
roll  vortices  must  appear  as  part  of  the  turbulence.  To 
examine  the  relative  contributions  of  the  smaller  scale  random 
turbulence  and  the  relatively  organized  large  scale  rolls,  a 
two  dimensional  calculation  was  performed  in  the  (y,z)  plane 
perpendicular  to  the  axis  of  the  rolls  (Leveller,  et  al.  1979). 
Periodic  boundary  conditions  were  applied  so  that  the  domain 
covers  a  wavelength  X  in  the  v  direction.  The  vertical 
domain  extended  from  the  surface  (represented  by  an  effective 
roughness  zQ)  to  slightly  above  the  capping  inversion.  The 
coordinate  system  permits  the  horizontal  roll  vortices  to  appear 
as  part  of  the  ensemble  mean  motion  when  calculating  the 
unsteady  flow  in  the  unstable  boundary  layer. 

The  numerical  code  partitions  the  energy  between  the 
mean  background  motion  which  is  a  function  of  the  vertical 
coordinate  only;  the  mean  ouasi-periodic ,  two-dimens ional ,  roll 
vortex  motion  which  is  a  function  of  z  and  y  ;  and  the  more 
random  turbulent  motion  which,  although  three-dimensional  in 
character,  is  only  a  function  of  y  and  z  in  the  mean. 

The  energy  in  the  organized  roll  motion  varies  with  the  ratio 
of  the  wavelength,  X  ,  of  the  roll  to  the  inversion  height 
z^  ,  the  instability  of  the  layer  as  measured  by  the  ratio  of 
the  Konin-Obukhov  length,  L  ,  to  z.  ,  and  the  angle,  a  ,  between 
the  roll  axis  and  the  geostrophic  wind.  Figure  10  illustrates 
the  cross-sectional  structure  of  the  stream  function  of  the 
roll  motion  for  an  angle  of  -10°  ,  X/z^  =  3  and  L/z^  =  -0.1 
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to  approximate  observations.  Surface  conditions  are  chosen 
appropriate  for  a  sea  surface  at  constant  temperature.  Results 
are  scaled  in  terms  of  the  characteristic  velocity 


to  minimize  the  influence  of  changes  in  w6  .  (The  angled 

s 

bracket  represents  an  average  over  y.)  The  geostrophic 
velocity  is  held  at  Uq  *  10  m/s  throughout  the  calculations. 

The  corresponding  cross-sectional  structure  of  the 
turbulent  kinetic  energy  and  the  humidity  are  shown  in 
Figs.  11  and  12.  The  numerical  results  for  the  relative 
contributions  of  the  mean  roll  motion  and  the  roll  nodulated 
turbulence  to  the  transport  of  momentum  and  heat  compare 
reasonably  (Lewellen  et  al .  iy79)  with  the  observations  of 
LeMone  (1976). 

In  Figs.  11  and  12  the  background  levels  have  been 
o 

removed;  hence  Aq  shows  some  regions  of  negative  values. 

The  updraft  region  shown  in  Fig.  10  produces  a  band  with  a 

spread  of  about  0.25  X  which  has  greater  than  background 
2  2 

C  and  q  .  The  rest  of  the  q  region  is  dominated  by 

2 

a  much  broader  area  of  less  than  average  values  of  Aq 
Figure  12  shows  a  noticeable,  almost  jet-like  character  to 
AC  .  The  heat  flux  is  concentrated  upward  in  the  same  region 
also,  producing  a  temperature  overshoot. 

The  vertical  structure  of  the  flux  of  a  species  such  as 
humidity  from  the  surface  is  shown  in  Fig.  13.  At  the  surface 
the  vertical  transport  is  all  carried  by  the  small  scale 
turbulence  designated  by  <wc>  but  at  z/z ^  :  0.4  approximately 

2/3  of  the  vertical  flux  is  contained  in  the  mean  roll  motion 
designated  by  <WC>  .  The  heat  flux  structure  is  similar 
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(Lewellen,  et  al.  1979}  except  for  an  undershoot  in  <W0> 
above  z/z ^  -  0.75  forced  by  the  capping  inversion. 

As  illustrated  in  Fig.  9,  the  most  common  way  of 
visualizing  the  organization  of  the  mixed  layer  into  a 
longitudinal  roll-vortex  structure  is  the  observation  of 
cloud  streets.  Figure  14  shows  the  relationship  between 
liquid-water-content  contours  and  streamlines  of  the  roll 
motion.  The  case  shown  in  Figure  14  is  somewhat  different 
from  that  detailed  in  Figs.  10-13.  In  the  calculation  for  the 
earlier  figures  the  species  C  transported  from  the  surface 
is  treated  as  a  passive  tracer.  For  the  calculation  of  Fig.  14, 
the  species  was  interpreted  as  humidity  and  allowed  to  condense 
under  appropriate  conditions.  The  subsequent  energy  release  is 
responsible  for  the  minor  change  in  streamline  appearance 
between  Fig.  10  and  Fig.  14. 

When  radiative  cooling  was  added  to  the  computation 
illustrated  in  Fig.  14,  the  cloudiness  increased  to  a  complete 
stratus  condition.  The  subsequent  strong  increase  in  cooling 
from  the  top  of  the  mixed  layer  forced  a  marked  change  in  the 
roll  structure.  The  wavelength  decreased  sharply  and  the 
resulting  mismatch  between  the  computational  domain  periodic 
conditions  and  the  desired  wavelength  of  the  phenomenon 
destroyed  most  of  the  organization  of  the  simulation. 


REFRACTIVE  INDEX  FLUCTUATIONS 


The  basic  turbulent  correlations  predicted  by  the  model 

can  be  used  to  predict  fluctuations  in  the  refractive  index 

for  the  propagation  of  electromagnetic  waves.  A  parameter 

of  major  concern  to  the  engineer  is  the  Refractive  Index 

2 

Structure  Function,  Cn  ,  which  describes  the  intensity  of 

refractive  index  fluctuations  over  the  inertial  subrange  of 

turbulent  eddies.  The  combination  of  primary  variables, 

temperature  and  humidity  variances  and  covariances  required 
2 

to  compute  Cn  depends  upon  the  wavelength  being  refracted. 
A  formula  for  radar  C?  was  derived  by  Lewellen  and  Teske 
(1975). 


where  B  =  7,730  ,  p  is  in  atmospheres,  T  is  °K  and  H 

is  water  vapor  mixing  ratio.  It  includes  the  influence  of 

mean  humidity  level  which  is  often  important  in  the  marine 

boundary  layer,  but  was  neglected  by  Wesely  (1976).  Burk 

(1979)  deals  with  the  prediction  of  acoustic,  optical  and 
2 

microwave  CR  throughout  the  planetary  boundary  layer 
using  a  model  (Burk,  1977)  similar  to  that  presented  here. 

2 

Figure  15  gives  the  contours  of  Cn  predicted  by  our 
model  for  one  set  of  conditions.  (Sea  surface  temperature  = 
20°C,  Ug  *  10  m/sec  ,  a  *  10  ,  z^/L  =  10,  30/9z  (above 
mixed  layer)  *  ,01°/m,  -  2.5  g/Kg.)  This  set  of  conditions 

differs  from  that  given  in  Figure  14  because  we  wished  to  show 
a  cloud-free  case.  The  figures  show  that  the  highest  values 
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of  C  occur  near  the  surface  and  in  the  vicinity  of  the 
n 

I  inversion.  There  is  also  considerable  variation  across  the 

roll. 

A  comparison  between  one-dimensional  and  two-dimensional 
temperature  and  humidity  variances  shows  that  the  major 
*  difference  between  one-dimensional  and  two-dimensional  pre¬ 

dictions  occurs  in  the  region  of  the  inversion. 
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CONCLUDING  REMARKS 


Viable  means  are  available  for  estimating  the  distribu¬ 
tion  of  water  vapor  and  its  variance  in  the  atmosphere  boundary 
layer.  However,  the  techniques  discussed  in  previous  papers 
in  this  volume  for  observing  water  vapor  in  the  atmosphere 
should  provide  the  means  for  obtaining  data  which  can  be  used 
to  improve  these  models.  This  is  particularly/  true  in  the 
neighborhood  of  the  capping  inversion  and  for  the  stable 
boundary  layer.  Although  a  number  of  simple  entrainment  models 
have  been  published  for  mixing  across  the  inversion  (e.g., 
Deardorff  1979  and  Yamada  and  Berman  1979)  none  of  these  are 
very  accurate  for  predicting  the  temperature  and  humidity 
fluctuations  which  occur  there.  A  need  exists  for  sufficiently 
complete  data  sets  which  can  be  used  to  either  provide  empirical 
correlations  of  mixing  across  this  layer  or  to  test  the  validity 
of  second-order  closure  models  such  as  the  one  described  herein. 

Comparisons  between  the  results  of  two-dimensional  and 
one-dimensional  simulations  of  the  same  mixed  layer  problem 
(Lewellen  et  al.,  1979)  lead  to  ideas  as  to  how  the  one¬ 
dimensional  model  may  be  improved  by  incorporating  some 
features  of  an  anisotropic  length  scale.  We  are  currently 
attempting  to  complete  this  model  modif ication .  Alternative 
approaches  to  improving  the  model  are  being  pursued  by  others 
(e.g.,  Lumley  et  al,  1978;  Wyngaard,  1979;  Gibson  and  Launder, 
1978) . 
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Abstract 

Recent  applications  of  the  turbulent  transport  model  originally  deseloped  by  Donaldson  to 
some  problems  of  practical  importance  in  micrometer  rology  are  discussed  Four  particular 
examples  considered  are  the  local  boundary  lax et  gust  front  emanating  ftom  a  thunderstorm;, 
the  low  -level  wind  and  turbulent  distributions  of  a  tornado,  the  transport  of  momentum  ,  heat 
and  species  withir.  a  surface  la>er  canopy ,  and  long;rudma!  roll  vortices  in  the  unstable  plane- 
tan  boundary  layer  Results  for  the  last  example  are  discussed  in  some  deiail  Comparison1 
are  made  between  a  one-dimenaona!  and  a  two-dimensional  computation  of  this  phenomenon 


List  of  Symbols 

C  Mean  species  concentration 

Species  fluctuation 
g  Gravitational  acceleration 

H  Height  of  vegetation  canopy 

h  Heighi  of  numerical  gnd 

L  Morun-Obukhov  length 

P  Mean  pressure 

q  Root-mean  velocity  fluctuation 

r  Radius 

f  Radius  at  which  maximum  swirl  occurs 

t  Tune 

Mean  velocity  components 
u,,u,.uK  Fluctuating  velocity  com¬ 
ponents 

Vc  Geostrophic  wind 

U,  Friction  velocity 

V  Mean  longitudinal  velocity 

u  Fluctuating  longitudinal  velocity 

V  Mean  transverse  velocity 

v  Fluctuating  transverse  velocity 

W  Mean  Vertical  velocity 

w  Fluctuating  vertical  velocity 

w.  Characteristic  velocity  for  free  con¬ 
vection 


x,  y  Horizontal  coordinates 
I  Vertical  coordinate 

r,  Height  of  mixed  lay  r 

:0  Surface  roughness  l...gh: 
q  Angle  between  the  roll  axis  and  the 

geostrophic  wind 
6,;  Kronecke r  delta 

e(/-i  Alternating  tensor 

A  Turbulent  macroscale 

X  Wavelength  of  the  roll 

©  Mean  temperature 

8  Temperature  fluctuation 

6.  -wS/T. 

ou  Standard  deviation  of  longitudinal 

wind  velocity 

oH  Standard  deviation  of  vertical 
wind  velocity 
v  Kinematic  viscosity 

p  Density 

^  Stream  function 

n  Earth's  rotation 
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Introduction 


The  basic  rurbulent  transport  model  used  in  these  applications  is  that  outlined  bv  DonsU son 
(1|  and  described  in  detail  bv  Lewetkn  |dj  Simple  closure  appropriations  were  made  to  the 
equations  for  the  ensemble-averaged  single-point  second-order  moments  of  the  fluctuating 
variables  The  theme  of  the  present  papet  is  to  discuss  the  various  applications  rather  than  to 
rede  rive  the  mode!  Since  invanance  is  an  essential  feature  of  the  model,  the  same  modeling 
coefficients  are  used  for  all  the  app’. cations  The  basic  modeled  equations  are  gnen  ir.  the  Ap¬ 
pendix  Although  this  simple  closure  cannot  be  expected  to  faithfullv  represent  alu  of  the  ir.for 
mation  present  m  complex  micrometeoro'ogical  turbulent  flows,  the  added  ph>  sics  incorporat¬ 
ed  ir  the  se.ond-order  correlation  equations  permits  an  influence  of  such  phenomena  as  buos  • 
ar.;;,  ard  rotation  to  he  calculated  directlv  without  resorting  to  ad  hoc  e dd\  viscosm  fixes 

The  H  rst  three  examples  w.ll  onlv  be  discussed  verx  bnefh  ,  but  the  fourth  will  be  pre¬ 
sumed  u  greater  detail  because  it  not  on!;-,  helps  to  illuminate  the  influence  that  roll  vortices 
eve-:  or  the  turbulent  trar  sport  of  momentum,  he  a1  and  species,  but  also  the  strong  ard  weak 
po.nts  of  the  basic  turbulence  closure  model 


Thunderstorm  Gust  Front 

TT  e  gust  front  caused  bs  the  cold  outflow  emanating  from  a  thunde'storrr.  ha-  beer,  modeled 
w.-j.  ar  axis.- tnmetric  implementation  of  ou’  turbulence  mode!  f.'j  The  sin  u. ’a  red  flow  field  is 
ti’cst'ated  in  Fig  1  \ke  ideaVe  the  outflow  from,  the  thunderstorm  a-  a  co'd  jet  of  tempera- 

tv  n  impir-ging  normal  to  the  ground,  released  at  a  height  -*m4>  with  vertical  velocin  kin 
T'  .•  ’-cmperarure  defect  be|0w  the  ambient  temperarure  is  caused  bv  evaporation  ol  falling  rain 
h>  relaroe.'j  dr\  air  at  some  altitude  ; 

The  primarx  sanables  in  this  simulauor.  are  the  temperature  defect  ( below  ambienii  of  die 
jet  it-  d.ameter.  and  the  heigh:  a’  which  it  is  re'ea-ed  The  larger  scale  updraft  wither,  wbu.-h 
the  downdraft  is  embedded  will  plav  a  role  in  retard. r.g  the  latetune  gust  front  within  die  com¬ 
putational  djmain.  but  should  not  be  a  crinca  facto;  a-  long  as  the  inflow  velocitv  is  smaller 
than  the  resu'ting  internal  simulated  velocities  The  other  critical  parameie r.  ,:'e  the  surfa.e 
temperature  the  surface  roughness  and  the  stahibtx  of  the  amb.ent  atmosphere  Ir.  an  effort 
tfi  esa'uate  ha.  ardous  conditions  for  ar.raf;  operauons.  a  sensilivirv  anaJv  s:s  ha-  been  per¬ 
formed  or.  the  gust  front  as  a  function  of  five  different  dimensionless  ph  x  sic  a;  parameters 

lJ-  M 


*  -  -  8  -  'Z’Z 


Fig  I.  Coordinate  <>  viem  for 
the  axi's  mmcinc  radial 
thundemonn  fu-t  front 
simulation 
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Fig  2.  KopJ'Th'C*  wOn'tur.i  temperature  defect  ' <?*  =  -  2  C*  ai  seural  time*  attct  ir.i!;jl;?aiipn 


F«  3.  Model  simulation  of  the  leading  edge  of  a  guM  front  *hen  visualized  in  term'-  of  temperature 
d.  fVc?  intenv»r> 


One  way  to  folio*  the  siructure  of  the  developing  gust  is  to  observe  the  movement  of  a 
fixed  temperature  value  line,  in  this  case  ©  =  -  2  °C,  as  shown  in  Fig  2  for  the  simulation 
sketched  in  Fig  1  At  t  -  0  the  ininal  linear  profile  is  shown,  but  by  t  *  200  s  the  structure 
has  developed  into  a  moving  front.  At  r  -  340  s  the  front  has  tom  away  fiom  the  downdraft 
region.  At  later  times  the  strength  of  the  gust  decreases  as  the  effect  of  the  area  change  be¬ 
comes  more  important.  However,  the  height  of  the  gust  appears  to  grow  slowly,  so  that  by 
t  •  1 ,000  s  the  ©  *  -  2  °C  line  reaches  neatlv  1  km  m  alutude  at  a  distance  of  nearly  8  km 
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Fij  4.  The  Ie3dmj  edge  of  a  thunder¬ 
storm  gusi  from  a-  rr.rWed  b\  d-v. 
(Photo  b>  Andre*  W anon. 

Florida  1975) 

from  the  source  centerline.  The  front  is  still  quite  sttong  it  this  point,  and  is  moving  outlaid 
at  approximately  2.6  m's. 

When  our  numerically  simulated  fiont  is  visualized  in  terms  of  the  temperature  defect 
intensity .  the  leading  edge  appears  as  shown  m  Fig  3.  There  is  a  very  strong  qualitative  sin-- 
Ian  tv  between  this  picture  and  the  gust  ftom  as  visualized  b>  dust  in  Fig  4  (also  see  picture 
in  (6]l  The  quantitative  predictions  of  maximum  mean  velocity .  velocity  shear  and  variance 
also  appear  consistent  wish  available  observations. 


Tornado  Boundary  Laver 

The  strongest  winds  occurring  in  nature  are  those  in  a  lornado.  in  order  to  estimate  wind  loads 
on  structures,  we  have  modeled  the  tornado  boundary  layer  [7]  again  using  tht  axisymmetnc 
version  of  our  basic  turbulent  transport  model.  Here  we  have  chosen  a  computer  simulation 
which  show  s  a  resemblance  to  tht  1974  Xenia.  Ohio  tornado  The  domain's  outer  radius  and 
top  are  both  placed  at  400  m.  The  surface  roughness  is  taken  to  be  large  (0  4  m )  to  correspond 
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Ff  5.  Moan  tangential  veloc.iv  contour-  *v  Fg  6.  Mean  radial  vetocit)  contour-  for  thr  same 

predicted  b>  an  a\is>  mmnric  tornado  mode!  conditions  a-  Fig  5 

Con'our»  labeled  in  tenth-  of  mavimum  f,  and 

radial  and  vertical  coordinate  normalutd  bs  the 

rad.u-  at  which  Vmil  occutv 
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to  a  relatively  urban  area  All  components  of  the  velocity  are  specified  at  the  outer  radius  in¬ 
flow  boundary  and  relatively  free  zero  slope  conditions  applied  at  the  top  of  the  domain. 

The  mean  tangential  and  radial  velocity  field  is  given  in  Figs  S  and  6.  The  velocities  are 
normalized  by  the  maximum  tangential  velocity  which  occurs  at  r  *  180  m.z  *  100  m.  Cur¬ 
rently  .  boundary  conditions  are  not  known  precisely  enough  for  the  absolute  value  of  the 
maximum  velocity  to  be  predicted  by  the  model  However,  from  photogrammetric  analy  ses, 
n  appears  the  tangential  velocity  should  be  scaled  by  a  maximum  value  of  approximately 
10C' m  s 

A  dominant  feature  of  the  flowfields  computed  to  date  (8, 9)  is  that  the  maximum  winds 
occur  at  quite  low  altitudes  The  surface  layer  plays  the  rather  paradoxical  role  of  increasing 
the  velocity  in  the  neighborhood  of  the  ground  at  small  radii,  before  actually  reducing  the 
velocity  to  zero  at  the  surface.  The  strong  radial  inflow  m  the  surface  lay  er  permits  the  stream¬ 
lines  to  penetrate  io  smaller  radii  here  than  at  higher  altitudes  Even  without  total  conserva¬ 
tion  of  angular  momentum,  this  permits  higher  swirling  velocities  to  develop.  The  simulation 
also  predicts  thai  the  rm  .s  average  of  the  velocity  fluctuations  reaches  values  as  high  as 
0.3  l'm»*  This  large  magnitude  of  the  fluctuating  velocities  indicates  that  they  should  be  con¬ 
sidered  m  setting  design  criteria  The  maximum  damage  will  most  likely  occur  where  the  fluc¬ 
tuating  velocities  add  to  the  mean  velocity  . 


Flow  Within  a  Surface  Canopy 

A  canopy  of  vegetation  presents  a  complex  lower  boundary  for  aimosphenc  flow  s  For  flow 
well  above  th.s  canopy .  n  is  usually  adequate  to  characterize  the  boundary  in  terms  of  only  an 
aerody  namic  roughness.  :0  But  for  flow  within  the  canopy  itself,  a  more  detailed  representa¬ 
tion  is  required  Our  motivation  for  detailing  this  flow  is  to  aid  in  the  prediction  of  the  dry 
deposition  of  gaseous  SO:  and  paniculate  sulfate 


7.  Comparison  of  A  R  A  P  model  prediction  with  data  from  Shffw  ft  a).  1 11,  12)  in  and  above  a  corn 
canop  where  H  is  canop'  height.  l\  the  friction  velocir  above  the  canop> .  V/V»  the  normalized  mean 
longitudinal  vclocin  ou'i\  and  ow/lr.  the  dimensionless  standard  deviations  of  the  longitudinal 
and  vertical  velocities,  and  t  -ww /L'J)  the  dimensionless  Res  nolds  stress 
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A  second-order  closure  model  for  canopy  flow  has  recently  been  presented  by  Ml  ton  and 
Shaw  [10].  The  principal  difference  between  their  model  and  ours  is  that  we  consider  heat 
and  species  transport  as  well  as  momentum  transport.  We  also  use  a  somewhat  more  general 
representation  of  the  drag  per  unit  -,'lume  of  the  vegetation.  The  comparison  of  the  mean 
distributions  of  velocity.  Reynolds  stress  and  velocity  variance  predicted  by  our  model  are 
compared  with  some  limited  data  for  a  com  canopy  [11,  12]  in  Fig.  7.  Although  detailed  plant 
area  densities  must  be  given  for  any  ty  pe  of  vegetation,  the  model  can  predict  the  variation  in 
surface  layer  heat  and  species  transport  as  a  function  of  surface  Reynolds  number  ( l',z0!v ), 
Prandtl  number  and  Schmidt  number  for  any  given  plant  area  density  distribution  [13]. 


Roll  Vortices  in  the  Planetary  Boundary  Layer 

The  final  example,  which  we  will  discuss  in  more  detail,  is  that  of  computing  the  longitudinal 
roll  vortices  which  often  occur  in  the  unstable  planetary  boundary  layer  [  14],  This  boundary- 
layer  feature  can  often  be  visualized  as  in  Fig.  8  by  the  streets  of  clouds  which  occur  at  the  top 
of  the  up  flow  region  between  the  roll  vortices.  The  simulation  is  performed  on  a  two-dimen¬ 
sional  grid,  stretching  over  a  wavelength  X  in  they  direction,  and  from  the  surface  (represented 
by  an  effecuve  by  drody  namic  roughness  z0)  to  the  total  height  h  of  the  numerical  grid  At  a 
position  Zj<h.  the  temperature  profile  is  damped  back  to  its  background  gradient,  capping 
the  domain  and  forcing  to  become  the  inversion  height.  All  mean  variables  and  turbulence 
quantities  are  initialized  with  profiles  from  a  compauble  boundary  layer  solution  of  our  one¬ 
dimensional  code.  The  spatially  homogeneous  one -dimensional  code  necessarily  treats  the  mo¬ 
tion  arising  from  the  roll  vortices  as  a  part  of  the  turbulence.  Periodic  boundary  conditions  are 
applied  to  the  two-dimensional  Cartesian  version  of  our  model.  This  permits  horizontal  roll 
vortices  to  appear  as  part  of  the  ensemble  mean  motion  when  calculating  the  unsteady  flow  in 
the  unstable  boundary  layer. 

The  numerical  code  is  free  to  partition  the  energy  between  the  mean  background  motion 
which  is  a  function  of  the  vertical  coordinate  only;  the  mean  quasi-penodic.  two-dimensional, 
roll  vortex  motion  which  is  a  function  of  :  and y: and  the  more  random  turbulent  motion  which, 
although  three-dimensional  in  character,  is  only  a  funcuon  ofy  and  z  in  the  mean  The  energy 
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F«.  8.  View  of  cloud  streets 
ovet  the  coast  of  Georgia  as 
seen  by  Apollo  1 1 4 ) 


Fig,  9.  Roll  stream  function 
contour-  for  the  example 
case  of  a  -  -IP.  \i:,  *  3. 
Ama» ■'-<  '  0. 1 .  The  profile- 
are  normalized  h>  a  maxi¬ 
mum  value  A.mjl's».r, 

*  0.28'  i  A.  =  v  -  -v  ’I. 
with  r9  denotinp  :90"  of 
maximum  value.  : "’O" . 
etc  The  arrow  -  indicate 
flow  direction 


in  the  organized  roll  motion  varies  with  the  ratio  of  the  wavelength.  X.  of  the  roll  to  the  inver¬ 
sion  height  r,.  the  instability  of  the  lay  er  as  measured  by  the  ratio  of  the  Momn-Obakhos 
length. L.  to  and  the  angle  between  the  roll  axis  and  the  geosnophic  wind  Figure  9  illus¬ 
trates  the  tx  pical  cross-sectional  structure  of  the  stream  function  of  the  roll  motion  for  an 
ai  cle  of  -  10'  .X  *  3  and  L  =  -0. 1.  Surface  conditions  are  chosen  appropriate  for  a  sea 
s-n'ace  at  constant  temperature.  Results  are  scaled  in  terms  of  the  characteristic  velocity 


i  ,,  ■  - u  -<  ,■ 

uo 

to  minimize  ihe  influence  of  changes  in  wf0 .  (The  angled  bracket  represents  an  average  over  v.) 
The  geostropb.ic  velocity  is  held  at  f  o  *  10  m  s  throughout  the  calculations.  To  maintain  the 
roll  structure,  we  unpose  periodic  boundary  conditions  on  all  variables  ai  v  =  ±X  2.  so  thai 
whaiever  goes  out  one  side  will  enier  the  other. 

The  corresponding  cross-sectional  structure  of  temperature,  velocity.  turbulent  kinetic 
energy .  and  passive  species  released  from  the  surface  perturbations  are  shown  in  Figs.  10 
through  13.  The  numerical  results  for  the  relative  contributions  of  the  mean  roll  morion  and 
the  toll  modulated  turbulence  to  the  transport  of  momentum  and  heat  compaje  reasonably  (15) 
with  the  observations  of  Lc.Uonc  [  16] 

Figures  10  through  13  show  strongly  similar  results  for  differences  above  background 
average  for  V.  ©.  q‘ ,  and  C.  Here  th*  background  levels  have  been  removed,  hence  A<?*  shows 
some  regions  of  negative  values.  What  is  obvious  from  all  four  plots  is  that  the  updraft  region 
shown  in  Fig  9  produces  a  band  of  less  than  background  6' and  greater  than  background  ©.  q'‘ 
and  C.  The  ry  pica!  spread  of  this  region  is  about  0.25  X.  The  rest  of  the  q1  region  is  dominated 
by  a  much  broader  area  of  less  than  average  values  of  A?2 .  Figure  7  shows  a  noticeable,  almost 
jet-like  character  to  V.  although  it  must  be  realized  that  Aly,,,,  =  0.096  <L9m,,  and  may  not 
be  that  easily  detected  Consistent  with  this  undershoot  in  A  V  is  an  overshoot  in  A©  (Fig  10). 
although  there  is  also  an  overshoot  coming  down  from  near  the  inversion  height  on  the  down¬ 
ward  side  of  the  roll.  A  very  broad  region  exists  where  A©  is  nearly  zero,  implying  little  varia¬ 
tion  of  0  iny.  These  results  are  all  consistent  with  LeMone's  observation  that  the  turbulence 
was  largest  where  the  vertical  flux  was  upward.  The  heat  flux  is  concentrated  upward  in  the 
same  region  also, producing  the  temperature  overshoot 


119 


Fig  13  shows  the  com  out  pattern  of  the  roll-averaged  AC,  showing  it  to  be  qualitative!) 
similar  to  the  transport  of  A©  off  the  surface.  No  return  transport  exists,  since  C  is  assumed 
to  go  to  a  small  value  at  the  top.  LeMone  observed  clouds  near  the  top  of  the  updraft  regions 
of  flow,  but  this  observation  should  not  be  confused  with  our  species  results.  The  passive  spe¬ 
cies  is  fairly  uniform  across  the  domain  ( ACmiX  =  0.050  <Omax ). 

If  the  closure  modeling  w’ere  exact,  the  one-dimensional  (horizontally  homogeneous)  anal¬ 
ysis  of  this  problem,  which  averages  over  the  periodic  large  eddy  roll  structure  and  treats  the 
roll  energy  as  part  of  the  turbulent  kinetic  energy  .  would  yield  the  same  answer  as  obtained  by 
horizontally  averagng  over  the  two-dimensional  result  Since  the  two-dimensional  computation 
allows  the  dominant  large  eddy  structure  to  be  determined,  the  closure  modeling  should  be  less 
critical  in  this  computation  than  it  is  in  the  analogous  one -dimensional  computation.  For  this 
phenomenon,  which  has  a  strong  two-dimensional  character,  out  mode!  represents  an  inter- 
medute  step  between  depending  completely  on  closure  modeling  and  depending  on  sub-grid 


!4.  Comparison  between  the  one -dim?*' mot, a  3^d 
two-dimensional  predictions  for  vertical  velocity  vari¬ 
ance.  For  the  ixso-dimentior  a!  calculator,  the  v- 
averaeed.  large  scaJc  roll  contribution  is  denoted  b\ 
<H<2-  and  the  smaller  scale  turbulence  ku 


Fig.  15.  Comparison  between  the  one-dimen^onal  and 
tv  o-dimenMonal  prediction*-  foj  tran^ci'C  velocity 
variance.  For  the  two-dimensional  calculation,  the  v- 
averaged.  large-scale  roll  contribution  iv  denoted  b\ 

<  V*  and  the  smaller  scale  turbulence  b\  vi' 


Fig.  16.  Companion  between  the  onc*dimensionj!  and 
two-dimensional  prediction^  for  longitudinal  velocit) 
variance.  For  the  tu c-dimensiona)  calculation,  the  y * 
averaged,  large  scale  roll  contribution  is  denoted  b\ 
<(/*  and  the  smaller  scale  turbulence  b>  iuu' 


closure  modeling  ( I  7j.  Thus,  the  two-dimensional  result  can  be  used  to  check  strengths  and 
weaknesses  of  the  one -dimensional  model. 

The  separate  components  of  the  velocity  variance  are  shown  in  Figs.  14  through  16.  There 
is  a  reasonable  comparison  between  the  vertical  variance  predicted  by  the  one-dimensional 
model  and  the  sp3tia!  average  of  the  two-dimensional  results.  But  the  other  two  components 
show  significant  differences. 
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The  analogous  comparisons  of  the  vertical  momentum  flux  and  the  vertical  heat  flux  are 
shown  in  Figs  17  and  18.  In  the  middle  of  the  mixed  lay er,  a  major  portion  of  the  turbulent 
transport  is  earned  by  the  large  roll  structure.  In  spite  of  this,  the  one-dimensional  model  is 
able  to  mode!  the  turbulent  transport  reasonably  well  except  in  the  vicinity  of  the  upper  inver¬ 
sion  Both  the  one-dimensional  and  two-dimensional  models  predict  aheat  flux  which  is  coun¬ 
ter  gradient  over  approximately  one-half  of  the  mixed  laser  depth.  The  two-dimensional  model 
predicts  a  larger  undershoot  in  the  heat  flux  near  the  inversion.  This  seems  to  be  associated 
with  the  fact  that  the  two-dimensional  model  can  account  for  more  influence  of  the  wave  mo¬ 
tion  in  the  stable  region  than,  is  represented  in  the  one-dimensional  model.  This  is  perhaps  shown 
better  in  Fig.  19  which  shows  the  temperature  variance.  By  far,  the  largest  contribution  to  tem¬ 
perature  variance  at  the  top  of  the  mixed  layer  is  the  wave  motions  induced  by  the  large  rolls. 
However,  even  the  small  scale  turbulent  variance  is  approximately  a  factor  of  2  larger  than  that 
predicted  by  the  one-dimensional  model  at  -  1 . 
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F^.  19.  Compandor  between  the  one-dimen  Mona]  and  ru  o-dimenwona)  prediction'  for  the  temperature 
variance  The ,) -averaged  roll  contribution  ^  denoted  b\  ©‘  and  the  turbulent  contribution  b> 


Figures  14  through  19  show  that  although  the  turbulent  transport  compares  quite  favorably 
for  these  two  models,  there  is  considerably  more  kinetic  energy  near  the  surface  and  near  the 
inversion  in  the  more  correct  two-dimensional  model.  The  deficiencies  of  the  one-dimensional 
model  appear  to  reflect  the  difficulty  of  modeling  with  an  isotropic  turbulent  macroscale.  It  is 
dear  from  the  two-dimensional  simulation  as  it  is  from  field  observations  [18]  that  at  both  the 
top  and  bottom  of  the  domain,  the  characteristic  horizontal  scale  of  the  turbulence  is  much 
larger  than  the  characteristic  vertical  scale.  It  appears  that  the  next  step  in  improving  the  model 
calls  for  incorporating  some  structural  shape  in  the  scale  representation  ]  )9). 
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Appendix 


Model  Equations 

The  model  equations  of  motion  fot  an  incompressible  fluid  in  the  presence  of  both  a  gravita¬ 
tional  and  a  Conolis  body  force,  with  the  mean  variables  denoted  by  capitals  and  the  turbulent 
fluctuations  by  lower-case,  may  be  written  in  general  tensor  notation  as  follows. 


5W,  “i  .  -  ^  -i  ^ 

8r  '  bx,  bx,  p  3x,  0O 


ail-  =  0 

3x, 


84  * 

or  ' 


30 

3x, 


bu,8 

3x, 


bu,u.  bu,u.  _  bl'.  _  bV,  u.8  ufi 

+  t*  =  -u,uk  —!■  -u,uk  - ~  +S,  -t  **>  -t 

bi  dx*  dx*  aX^  ©o  ©o 

-  -t/sA  U,U,  -  -  C;H 

3  du.ti,  o  i _  .  a‘  ,  ai 

,0.3  -  •«  A  -j—t  |  -  -  K«,  -6„  3 , 


axk 


bu,6  bii,P  -  30  — v  31',  0*  n  — « 

'at  *  V>  "5T  =  a-~  +«-  51  ' 


3x 


3x 


ax,  “  e0 


.  ,  3  .  .  3u,e  1  0.75q  — 

403  s,  l,A  K1"  "a' 


r,  lf:  .-:.7  ±  I, A  it,  - 

Of  dxj  bx ,  ox f  oxf  A 


£  ♦  l,  £  -  0.35  4 
dr  7  Ox,  qJ 


+  0.?5q  +  0.3  —  |?  A  3  ^ 


3x 


Ox, 


3x, 


_  °_?7J  +  °_8  .  Hi? 

q  bx,  I  q3  '  0Q 


The  overbar  represents  an  ensemble  average,  and  q3  =  u,u, . 

Along  with  the  differing  boundary  conditions  for  the  individual  problems,  an  upper  bound 
is  placed  on  A  depending  upon  the  spread  of  the  region  of  turbulence. 

For  the  canopy  flow  example,  variation  with  respect  to  one  coordinate  only  is  permitted. 
Also  extra  sink  terms  must  be  added  to  the  equations  for  the  mean  variables,  and  both  source 
and  sink  terms  added  to  the  Rey  nolds  stress  and  heat  Tux  equations  (13).  The  other  three  ex¬ 
amples  permit  dependence  of  variables  on  two  coordinates.  The  flow  may  be  eithet  planar  or 
axisymmetric. 
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